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ABSTRACT 

Transport phenomena in porous media are commonplace in our daily lives. Ex- 
amples and applications include heat and moisture transport in soils, baking and 
drying of food stuffs, curing of cement, and evaporation of fuels in wild fires. Of 
particular interest to this study are heat and moisture transport in unsaturated soils. 
Historically, mathematical models for these processes are derived by coupling clas- 
sical Darcy's, Fourier's, and Fick's laws with volume averaged conservation of mass 
and energy and empirically based source and sink terms. Recent experimental and 
mathematical research has proposed modifications and suggested limitations in these 
classical equations. The primary goal of this thesis is to derive a thermodynami- 
cally consistent system of equations for heat and moisture transport in terms of the 
chemical potential that addresses some of these limitations. The physical processes 
of interest are primarily diffusive in nature and, for that reason, we focus on using 
the macroscale chemical potential to build and simplify the models. The resulting 
coupled system of nonlinear partial differential equations is solved numerically and 
validated against the classical equations and against experimental data. It will be 
shown that under a mixture theoretic framework, the classical Richards' equation for 
saturation is supplemented with gradients in temperature, relative humidity, and the 
time rate of change of saturation. Furthermore, it will be shown that restating the 
water vapor diffusion equation in terms of chemical potential eliminates the necessity 
for an empirically based fitting parameter. 
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1. Introduction 

Water flow, water vapor diffusion, and heat transport within variably saturated 
soils (above the groundwater and below the soil surface) are important physical pro- 
cesses in evaporation studies, contaminant transport, and agriculture. The mathe- 
matical models governing these physical processes are typically combinations of clas- 
sical empirical law (e.g. Darcy's law) and volume averaged conservation laws. The 
resulting equations are valid in many situations, but recent experimental and math- 
ematical research has suggested modifications and corrections to these models. The 
primary goal of this thesis is to build a thermodynamically consistent mathematical 
model for heat and moisture transport that takes these recent advancements into con- 
sideration. To realize this goal Hybrid Mixture Theory (HMT) and the macroscale 
chemical potential are used as the primary modeling tools. 

1.1 Previous Work 

In 1856, Henri Darcy published his research on the use of sand filters to clean 
the water sources for the fountains in Dijon, France. He found that the flux of water 
across sand filters was directly proportional to the gradient of the pressure head. This 
simple observation has become known as Darcy's law and is one of the main modeling 
tools in hydrology and soil science [31J. Darcy's law is an example of a historical rule 
(or law) that has perpetuated to the present day. Darcy's law was originally derived 
for saturated porous media, but near the turn of the century it was extended for use 
in unsaturated soils. In 1931, L.A. Richards coupled Darcy's law with liquid mass 
balance to derive what is now known as Richards' equation. This equation relies on 
the assumption that Darcy's law is valid for unsaturated media, but it also relies on 
an empirically-derived relationship between pressure and saturation. 

The pressure-saturation relationship is known to be hysteretic in nature (de- 
pends on the direction of wetting), and only recently have researchers been able to 

move toward functional relationships that capture this effect [IE]. Correction terms 
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in Richards' equation have been proposed via HMT that suggest that the rate at 
which the capillary pressure is changing may play a role in the overall dynamics of 
the saturation [321 SZ] • This proposed, third-order, term in Richards' equation has 
only recently been studied mathematically and experimentally. One possible physical 
interpretation of this term is that is accounts for how fast the liquid-gas interfaces 
are rearranging at the pore scale. 

In the 1950s, Philip and deVries published their works on vapor and heat transport 
in porous media [321 El] • Their approach accounts for water flow in both the liquid 
and fluid phases in response to water content and temperature gradients in the soil. 
In order to account for the observation that Fickian diffusion inadequately describes 
diffusion in porous media, Philip and deVries implemented an enhancement factor, rj, 
to adjust the diffusion coefficient. This factor is fitted to the measured diffusion data 
for a particular medium. In [23] , Cass et al. found that rj increases with saturation 
(with r\ « 1 for dry soils). Philip and deVries proposed that thermal gradients and the 
condensation and evaporation through "liquid islands" are the pore-scale mechanisms 
that cause the observed enhancement. Counter intuitively, the governing equation 
holds at the macroscale while these mechanisms are inherently pore-scale. The Philip 
and de Vreis model has not been validated in a laboratory setting. 

More recent works question the validity of the Philip and deVries model. Shokri 
et al. [73] suggests that the coupling between the water flow and Fickian diffusion 
is the key to estimating the vapor flux. They suggest that under this consideration 
there is no need for the enhancement factor. Webb [80], and more recently Shahareeni 
et al. [70], showed that enhancement can exist in the absence of thermal gradients. 
It was initially thought that enhancement couldn't occurin the absence of thermal 
gradients and was therefore ignored. Based on the observtion that enhancement can 
occur without the need for thermal gradients, it is clear that the Philip and de Vries 
model needs modification. Cass [23] and Campbell [23] give a functional form of the 



enhancement factor that is commonly used (e.g. [§B]|78J), but relies on an empirical 
fitting parameter. In the present work we take the view of Shokri et al. that there is 
no need for the enhancement fact, and wederive a diffusion equation simply based on 
liquid and vapor flow. The novelty of the present approach is the use of the chemical 
potential as the driving force for both types of flow. 

For energy transport, the 1958 deVries model [32] is still commonly used (e.g. 
[68| 178]). Similar to the enhanced diffusion model, deVries built this model so as 
to account for the flux of the fluid phases. This is sensible as the fluid phases will 
certainly transport heat. More recently, Bennethum et al. [H] and Kleinfelter [52] 
used Hybrid Mixture Theory to derive heat transport equations in porous media 
(Bennethum et al. studies saturated porous media and Kleinfelter studied multiscale 
unsaturated media) . They verified many of the findings by deVries but also proposed 
several new terms associated with the physical processes of heat transport. In this 
work we extend the Bennethum et al. and Kleinfelter approachs to unsaturated 
media. 

1.2 Hybrid Mixture Theory and Thesis Goals 

To build the models in this work we make extensive use of Hybrid Mixture The- 
ory (HMT). HMT, statistical upscaling, and homogenization have all been used as 
techniques to re-derive, confirm, and extend Darcy's, Fick's, and Fourier's laws in 
porous media. For a technical summary of some of these methods see [30]. HMT is a 
term for the process of using volume averaged pore-scale conservation laws along with 
the second law of thermodynamics to give thermodynamically consistent constitutive 
equations in porous media. The technique as applied to porous media was developed 
by several parties, the most notable being Hassanizadeh and Gray J39J S3] and Cush- 
man et al. [TTI [T2"| |2"9"] . but the general principles were developed by Coleman and Noll 



In the present work we use HMT to derive new extensions to these laws in the 
case of unsaturated porous media. These extensions are then used to derive a model 
for total moisture transport in unsaturated soils. While this sort of modeling has 
been done in the past, rarely have the three principle physical process (movement 
of saturation fronts, vapor diffusion, and thermal conduction) been considered from 
first principles and put on the same theoretical footing (HMT in this case). No 
known work attempts to couple these three different effects together with one physical 
measurement: the chemical potential. This is one of the unique features of this work. 

In the most general sense, the chemical potential is a measure of the tendency of 
a substance (thinking particularly of a fluid or species) to diffuse. This diffusion could 
be of a species within a mixture (e.g. water vapor diffusing through air) or it could 
be a phase diffusing into another (e.g. water into a Darcy-type sand filter). The fact 
that most physical processes in porous media are of this type gives an inspiration for 
the potential usefulness of the chemical potential as a modeling tool. That is, from a 
broad point of view, it should be possible to restate the physical processes of moving 
saturation fronts and vapor diffusion more naturally by the chemical potential. This 
approach has not been thoroughly explored in the past since the chemical potential 
is not directly measureable and the theoretical footings of upscaling the chemical 
potential are relatively new. In the saturated case, extensions to Darcy's law have 
also been developed via HMT, and the results indicate that the macroscale chemical 
potential is a viable modeling tool for diffusive velocity in saturated porous media 
[T5l 1691 [81] . In the present work we give chemical potential forms of Darcy's and 
Fick's law as well as presenting simplifications to Fourier's law based on the chemical 
potential. In the case of a pure liquid phase we will show that the chemical potential 
form of Darcy's law is no different than the more traditional pressure formulation. In 
the gas phase, on the other hand, we will show that the pairings of chemical potential 
forms of Fick's and Darcy's laws gives a new form of the diffusion coefficient that 



does not need the enhancement factor indicated in the work by Phillip and DeVries 
[ST] . The chemical potential will finally be used to derive a novel form of Fourier's 
law for heat conduction in multiphase media. 

Once we have derived new forms of the classical constitutive equations we pair 
these equations with volume averaged conservation laws to give a coupled system of 
partial differential equations governing heat and moisture transport. The second law 
of thermodynamics is used to suggest additional closure conditions for each of the 
equations. Together, the system consists of a nonlinear pseudo-parabolic equation 
for saturation, a nonlinear parabolic equation for vapor diffusion, and a nonlinear 
parabolic-hyperbolic equation for heat transport. 

In summary, this work serves several purposes: (1) it is a step toward better 
understanding the role of the chemical potential in multiphase porous media, (2) it 
makes strides toward understanding the phenomenon of enhanced vapor diffusion in 
porous media, and (3) finally we propose a novel coupled system of equations for heat 
and moisture transport. 

1.3 Thesis Outline 

In Chapter [2] we take a step back from porous media and discuss pore-scale 
diffusion models. This is done in an attempt to elucidate the assumptions, derivations, 
and models used in various disciplines as there tends to be confusion about where the 
miriad of assumptions are valid. In this chapter we give mathematical and physical 
reasons for the many commonly used assumptions as well as proposing an alternative 
advection-diffusion model as compared to the popular Bird, Stewart, and Lightfoot 
model [18]. 

In Chapter [3] we present the necessary background information in order to un- 
derstand volume averaging and the exploitation of the entropy inequality. Much of 
this chapter is paraphrased from previous works, such as [UJ, [12], [39], 03], ED, ES]. In 

the beginning of Chapter [4] we use the tools from Chapter [3] to build and exploit a 
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version of the entropy inequality specific for multiphase media where each phase con- 
sists of multiple species. The remainder of Chapter [4] is dedicated to the exploitation 
of the entropy inequality for a novel choice of independent variables describing these 
media. Appendix [C] serves as a companion to this discussion as it gives the abstract 
formulation and logic of the entropy inequality. Throughout Chapter |4j the goal is to 
derive new forms of Darcy's, Fick's, and Fourier's laws and to propose extensions to 
these laws in terms of the macroscale chemical potential. As part of these derivations 
we arrive at new expressions for the pressure and wetting potentials in unsaturated 
media. All of these derivations are done in a general sense with as few assumptions 
as possible. This leaves open the possibilities of future research. 

In Chapter [5] we couple the results found from the exploitation of the entropy 
inequality (Chapter 111) with the volume averaged conservation laws derived in Chapter 
[3j In Section 5.1[ a more in-depth historical perspective of the classical equations used 



for heat and moisture transport is given to orient the reader to the recent research. 



Fluid transport equations are presented in Section 5.3.1 along with a discussion of 
the relationship between mass transfer and chemical potential. Considerable effort 
is put toward deriving a heat transport equation with the final equation presented 



in Section 5.3.2.3 In Section 5.4 several simplifying assumptions are presented in 



order to close the system of equations. In particular, Sections 5.4.1, 5.4.2, and 5.4.3 



give simplifications, assumptions, and dimensional analysis for the liquid, gas, and 



heat equations respectively. In Section 5.4.4 we present the remaining constitutive 



equations necessary to close the system of equations. Since so many assumptions and 
simplifications are made throughout the chapter a summary of all of the results is 



presented in Section 5.5 



In Chapter [6] we examine the proper regularity and assumptions needed for exis- 
tence and uniqueness of solutions. These results are preliminary and do not constitute 
a complete existence and uniqueness study for these equations. 
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In Chapter [7] we perform numerical analysis on the equations derived in Chapter [5] 



In Sections 7.1, 7.2, and 7.3 we examine numerical solutions and parameter sensitivity 



for the saturation equation, vapor diffusion equation, and the coupled saturation- 



vapor diffusion equations respectively. In Section [7T4| we compare numerical solutions 
to the fully coupled heat and moisture transport model to the experimental data 
collected in [78]. 

In Chapters [5] - [7] we work toward building and analyzing the saturation, vapor 
diffusion, and heat equations. The flow of thought for these chapters is to apply each 
set of new assumptions or simplifications to each of the three equations before moving 
to the next set of assumptions. That is, if a set of assumptions are proposed then the 
subsequent sections will apply those assumptions to the saturation, vapor diffusion, 
and heat equations in turn. Only then will the next set of assumptions be discussed. 
This is done so that each set of assumptions are only stated once and since many of 
the assumptions create interleaving effects between the equations. 

Finally, as an aid to the reader there are several appendices. Appendix |A~[ contains 
a nomenclature index for the pore-scale diffusion processes considered in Chapter [2] 



Appendix |B.1| contains a nomenclature index for the macroscale results in the re- 
maining chapters. There is some overlap between the nomenclature for these distinct 
parts, and effort has been made to not create any excessive notational confusions 



(even though this work is necessarily notation heavy). Appendix B.2 gives a list (in 
alphabetical order) of the upscaled definitions of variables defined in chapters [3] and |ij 
As mentioned previously, Appendix [C] gives an abstract view of the entropy inequal- 
ity in an effort to make the exploitation process more clear to the interested reader. 
Appendix [D] gives a summary of the results extracted from the entropy inequality 
in Chapter [4} This is done for ease of reference mostly on the author's part, but it 
is also done to provide an index of these results for use in future research. Finally, 
Appendix [E] gives several tables of dimensional quantities used throughout. 



It is suggested that the detail-oriented reader have Appendix [A] at hand when 
reading Chapter [2] and Appendix B.l at hand when reading chapters [3] - [5j There 
are some minor abuses of notation, but effort was made to bring them to the reader's 
attention whenever possible and to use notation that didn't confuse the immediate 
discussion. 



2. Fick's Law and Microscale Advection Diffusion Models 

This chapter consists of a short technical note related to pore-scale diffusion prob- 
lems. Vapor diffusion in macroscale porous media is an important phenomenon with 
many applications (e.g. evaporation from soils, moisture transport through filters, and 
C0 2 sequestration). In order to better understand macroscale diffusion it behooves 
the researcher to first understand pore-scale mechanics and models. This chapter at- 
templts to elucidate the models and assumptions used for diffusion at the pore-scale 
so that when we turn our attention to macroscale diffusion we are firmly grounded. 
A secondary goal of this chapter is to give a thorough discussion of the diffusion coef- 
ficient used in Fick's law. This is necessary since this coefficient is typically wrongly 
assumed constant for all choices of dependent variables (mass concentration, molar 
concentration, chemical potential, etc.). 

To make matters simpler, we focus our pore-scale discussion on the, so called, 
Stefan diffusion tube problem. This is a well-studied problem that models the dif- 
fusion of a species through an ideal binary gas mixture above a liquid-gas interface 
P HEJ 1221 E31 Ell EH ES]- This is an idealization of the juxtaposition of phases in a 
capillary tube geometry, and a capillary tube geometry is an idealization of geometry 
of pore-scale porous media. To derive a mathematical model for the time evolution 
of the evaporating (or condensing) species, one typically couples Fick's first law with 
the mass balance equation. 



In Section 2.1 we discuss the various forms of Fick's law and briefly discuss the 



relationships between the diffusion coefficients. In Section 2J2 we derive the tran- 
sient diffusion equations associated with Fick's law and compare with the associated 
equation of Bird, Stewart, and Lightfoot [18] (henceforth referred to as BSL). 

2.1 Comparison of Fick's Laws 

For a system consisting of an ideal mixture of water vapor, g v , and inert air, 

g a , Fick's law can be written in terms of molar concentration, mass concentration, 
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Table 2.1: Mass and molar flux forms of Fick's law 



Flux Type 


Flux Expression 


Fick's Law 


mass flux [ML- 2 T} 


J9j = p 9o v 9j,9 


J9 p i = -p9D p VC 9 * (2.1) 


molar flux [(mol) L~ 2 T] 


JBf = c^vl 3 ' 9 


JBi = -c 9 D c Vx 9 > (2.2) 


mass flux [ML~ 2 T] 


J9J = p 9j v 9j,9 


J %> = D ^(r 9jT ) V ^ (2 - 3) 


mole flux \(mol)L~ 2 T] 


J* c = dHyV* 


J%c = -D„,c (J^ VnV (2-4) 



or the chemical potential. This is potentially confusing since there are inherently 
different diffusion coefficients for the different forms of Fick's law. The purpose of 
this subsection is to clarify the relationships between these coefficients. In porous 
media it is common to use mass flux for Fick's law, but in chemistry (and related 
fields) it is more common to use molar flux. As such, we will make most of our 
comparisons between mass and molar flux. 

According to BSL |18j, the mass and molar forms of Fick's law are given by 
equations (O) and (O) (modified from BSL Table 17.8-2). In Table [O v 9 ^ 9 = 



v 9j _ v 9 j s th e diffusive velocity relative to a mass weighted velocity, Vc J ' 9 = v 9: > — v 9 
is the diffusive velocity relative to a mole weighted velocity, p 9i is the mass density 
of species j, C 9i = p 9j / p 9 is the mass concentration of species j in the mixture, 
c 9j = mol(gj)/vol(g) is the molar density of species j, and x 9j = c 9j /c 9 is the molar 
concentration of species j in the mixture. 

The chemical potential forms of Fick's law can be given in terms of two differ- 



ent types of chemical potential: mass weighted (equation (2.3)) or mole weighted 



(equation (2.4)). In physical chemistry and thermodynamics [211 ES] the chemical 
potential is known as the tendency for a species to diffuse, and for this reason it is a 
natural candidate for the statement of Fick's law (an exact thermodynamic definition 
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will be presented in subsequent chapters). In equations (2.3) and (2.4), ji 9 c 3 is the 
mole weighted chemical potential [J(moZ) _1 ] and \Tp is the mass weighted chemical 
potential [JM~ 1 ]. 

The reader should first note that the two fluxes are measured with respect to 
different velocities. The mass weighted velocity is p 9 v 9 = 'Yl,j =v a p 93 v 9: > an< ^ ^he m °l e 
weighted velocity is c 9 v 9 c = ^2j =va c 9j v 9:i where v 9i is the velocity of the species 
relative to a fixed coordinate system. This means that that there cannot be a direct 
comparison between the two different types of flux without considering them relative 
to the same frame of reference. Using these definitions of v 9 and v 9 we see that the 
difference between the two bulk velocities, v 9 — v 9 is 

JV 



V 9 



-)<:lj 



x 9j p& 






In (2.5), the summation over j indicates that this is an accumulation over the N 
species in the gas mixture. In future work we will be interested in the diffusion of 
water vapor (j = v) and will consider the gas mixture as binary: j e {v, a}, where 
j = a represents the mixture of all species that are not water vapor. Therefore we 



can write (2.5) as 



p9j _ x a 'ffi 



v 9 - v 9 = > 1 r — V 9i 



^ \ P 

j=v,a 

= (x 9a C 9v - x 9v C 9a ) v 9v + (x 9v C 9a - x 9a C 9 ")v 9a . (2.6) 

Converting to a mass weighted velocity we see that Vc J ' 9 = v 9j ' 9 + {v 9 — v 9 ), and 
therefore the molar flux is J 9 C 3 = c 9i v 9 c 3 ' 9 = c 9j v 9j ' 9 + c 9j (v 9 — v 9 ). Assuming that 
c 9j x 9k(j9i <^ i we see tha^ the difference between the frame of reference is potentially 
quite small. 

Next note that the diffusion coefficients are (initially) assumed to be different for 
each choice of independent variable as indicated by the subscripts. To compare D p 
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and D r we note that J 9 J = m 9j J 9 J and 



VC* = ( ^^ 2 ) Vx* (2.7) 

[x 9 im g i + x 9k m 9k ] 



to conclude that 



— J D P = D c , (2.8) 

where mP 3 is the molar mass of species j and the minuscule, k, represents the other 
species. If the molar density form of the diffusion coefficient, D c , is assumed to be 
constant (at constant temperature) we conclude that the mass density version of the 
diffusion coefficient is not constant (and visa versa). The fraction, C 9k /x 9k can be 
interpreted as the ratio of the molar mass of species k to the molar mass of the 
mixture. If j — v is the water vapor in an air-water mixture then k = a is the inert 
air and the scaling factor between the diffusion coefficients is the ratio of molar mass 
of the air to the molar mass of the mixture. For sufficiently dilute systems (where 
the amount of water vapor is small) the ratio is approximately 1 and the diffusion 
coefficients can be considered as approximately equal. 

For ideal air-water mixtures, the densities are related through p 9 = p 9v + p 9a and 
the water vapor density is related to the relative humidity through p 9v = p 9 s v at tp. Here 
we are taking p 9 / at as the saturated vapor density and tp as the relative humidity. At 
standard temperature and pressure we note that p 9 s v at ~ 0.02kg/m 3 and p 9 ps lkg/m 3 . 
This indicates that at standard temperature and pressure we can likely assume that 
the mixture is always sufficiently dilute. Therefore, in the systems under consideration 
we can assume that the diffusion coefficients are approximately equal. 

For the diffusion coefficients associated with the chemical potential forms of Fick's 
law we first observe that if we multiply and divide the right-hand side of the molar 
form by the molar mass of species j then 

*--M©p™"@?)G&K (2 - 9) 
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Here, R 9] = R/m 9] is the specific gas constant, and we have used p 9j = p 9 c 3 /m 9i . 
Again noting that J 9 * = rn 9] J 9 ^ c we conclude that D^ c = D^ p . 

It remains to compare D p to D pp and D c to D^ c . We focus here on the mass 
fluxes without loss of generality. If the mass fluxes are equal, then in particular 

Rearranging, it can be seen that 

D P VC* = C*-D„V [^- T 

(assuming constant temperature). From physical chemistry [56], recall that the chem- 
ical potential is related to a reference chemical potential (/if") and the ratio of partial 
pressure, p 9v , to bulk pressure, p 9 , via 

fjth> = /if + /^Tln ( ^7 ) • (2.10) 

Therefore, 

Using Dalton's law for ideal gases, p 9 = p 9v +p 9a , and using the specific gas constants 
we note that the partial pressure of species j can be written as p 93 = R 9j Tp 9i . 
Therefore, 

p 9 = R 9 -Tp 9v + R 9a Tp 9 % (2.12) 



and, after simplifying, 



f\9v s}9v 

(2.13) 



P 9 p 9 * + (S^)p 9a ' 

From the values found in Appendix py we see that R 9a /R 9v « 0.6 and therefore 



equation (2.13) is similar, but not equal to, the mass concentration, C 9v = p 9v /(p 9v + 



p 9a ). Defining C 9v = p 9v jp 9 we see that 



D p ~ \C 9 *) VCff-' { ' 
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where division is understood component wise (that is, equation (2.14) represents three 



equations when the gradient is understood in three spatial dimensions). The right- 



hand side of equation (2.14) is not constant at 1 for all densities, but the variation 
in the right-hand side depends mostly on the variation in p 9a in the gas mixture. 
Fortunately, the water vapor density is much smaller than the air-species density, 
and hence p 9a is approximately constant. In one spatial dimension, the right-hand 



side of (2.14) can therefore be approximated by 

/ x \ d I x \ 

r/ \ ._ \x+l) dx Vx+0.6/ 

a \ X ! -~ ( x \ _d_ ( x \ 

Vz+0.6/ dx \x+l) 

(where x = p 9v and p 9a ~ 1). It is easy to show that 0.98 < d(x) < 1 for < x < p g s v at . 
Furthermore, for sufficiently dilute mixtures, d(x) ~ 1, and we therefore conclude 
that D p « Dp. 

The conclusion from this subsection is that while the diffusion coefficients for the 
molar and mass flux forms of Fick's law are not the same, for dilute mixtures they 
can be approximated as equal. 

2.2 Transient Diffusion Models 

In porous media there is a phenomenon known as enhanced vapor diffusion }61j . 

This phenomenon states that vapor diffusion in porous media occurs faster than as 
predicted by Fickian diffusion models. This is merely a statement about the observed 
imbalance between Fickian diffusion and experimental measure. Since the ultimate 
goal of this work is to develop macroscale advection diffusion models, we seek to 
understand the pore-scale diffusion models so that in subsequent chapters we can 
tackle the enhanced diffusion problem. 

To build a transient model for molecular diffusion we couple Fick's law with the 
appropriate form of the mass balance equation. In the previous subsection we showed 
that the various forms of Fick's law are approximately equal (for sufficiently dilute 
mixtures), so the results stated here will only be in terms of the mass flux form of 



Fick's law (equation (2.1)). 
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The mass balance equation for species j in the gas phase can be written as 

dp 9 * 



dt 



+ V •{pPiv a ') = f 3 i (2.15) 



where v 9i is the velocity of species j within the gas mixture relative to a fixed frame of 
reference, and r 9i is a mass exchange term accounting for chemical reactions between 
species (HUES]. In the present work we assume that no chemical reactions occur, and 
therefore f 9j = 0. The combination of the mass balance equation with the mass flux 
form of Fick's law (for j = v) gives a transport equation for the mass of water vapor 
via advective, p 9v v 9 , and diffusive, J 9v , fluxes: 

^ + V -(J 9v +p 9 *v g ) = 0. (2.16) 

Substituting the mass flux form of Fick's lawp] (from equation (2.1)) we get 

^ + V-(p 9 V)=DV-(/) 9 V^). (2.17) 

Notice here that the diffusion coefficient has been factored out of the divergence 
operator. This is only valid in constant temperature environments. If the gas-phase 
density were constant in space then we would arrive at the traditional advection 



diffusion equation (by dividing equation (2.17) by p 9 of by rewriting the diffusion term 



as DV • Vp 911 ) and would need an expression for the bulk velocity in terms of density 
(or concentration) to close the equation. Unfortunately, if the density of the water 
vapor is allowed to vary then the density of the gas varies. Again, for sufficiently 
dilute mixtures the variation in gas-phase density is very small and the nonlinear 
diffusion on the right-hand side can be approximated by the linear diffusion term 
.DV ■ Vp 9 ". It should be noted here that this later case is what is typically thought 
of as "Fick's law" and is what leads to the traditional linear diffusion equation (when 
the advection term is neglected) 



lr The subscripts on the flux and the diffusion coefficient have been dropped since all of the versions 



presented in Table 2.1 are approximately equal 
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A different form of equation (2.17), suggested in BSL [IB], is derived by consid- 
ering the mass weighted bulk velocity. In a binary system, 

p 9v v 9v = p 9v v a„ 9 + p 9v v g = p g, v g v , 9 + p 9v (C 9 "v 9v + C 9a v 9a ) . 



Solving for p 9v v 9v 



p 9v v 9v 



P 



9v 



1 - o 



y9v,9 I n9v y9a 



(2.18) 



Using Fick's law for p9vv 9v ' 9 , and eliminating p9 v v 9v in (2.15) with (2.18) gives 

Q p 9v 



dt 



+ V ■ (p 9v v 9a ) = V 



Dp 9 
l-O 



VC fft 



(2.19) 



Whitaker [85] suggested that "one can develop convincing arguments in favor of ... " 
neglecting the air-species flux term. Certainly at steady state we can assume (as is 
done in BSL) that v 9a is approximately zero at the interface since "there is no net 
motion of [water vapor] away from the interface" [18j, but in the transient case this 
would constitute a change of frame of reference. This new frame of reference would 
be such that the inert air molecules are viewed as stationary with the water vapor 
diffusing through them. 



In either equation (2.17) or (2.19) one must find appropriate conditions or equa- 



tions to either neglect or rewrite the advective term, p 9v v 9 or pP v v 9a respectively. 
Typically this term is neglected in a pure diffusion problem. As these are two differ- 
ent simplifications of the same equation one must have different reasons for neglecting 
the advective term. The easiest fix for this issue is to couple with either the bulk 
gas mass balance equation or the air-species mass balance equation and to use the 
mass- weighted velocity: p 9 v 9 = Ylj=i p 9i v 9i . The point being that one cannot simply 
neglect the advection term in the transient case of either equation without proper con- 
sideration of the implications: a fixed bulk velocity or a changing frame of reference 
respectively. 



16 



A final comment can be make regarding equations (2.17) and (2.19). The bulk 
density term, p 9 , on the right-hand side of these equations is often factored out of the 
divergence operator. This is an error committed by several researchers [T8l [22l |5T] . 
The reasoning for assuming that the density is constant (and hence returning to a 
linear diffusion model in the absence of advection) is that in an ideal gas, p 9 = p 9 R 9 T. 
Under constant temperature conditions, and if the pressure is assumed constant, then 
the density is assumed be constant. There are two possible mistakes here. (1): If the 
species densities are allowed to vary then the bulk density must vary. (2): The value 
of R 9 will vary with the changing composition of the mixture (since the molar mass 
of the mixture changes). The effect of this is that, while the pressure may remain 
constant, the component parts are not necessarily constant and therefore cannot be 
factored from the divergence operator. 

2.3 Conclusion 

In this chapter we have compared various forms of Fick's law for molecular diffu- 
sion. We have shown that, while the diffusion coefficients are indeed different, under 
certain common circumstances the diffusion coefficients can be considered as approx- 
imately equal. It is common to take the diffusion coefficient as constant (or only a 
function of temperature), and in many cases it is safe to assume the same diffusion 
coefficient may be used in the common forms of Fick's law. 

In the transient case there are two natural formulations for (the mass flux form 
of) Fick's second law. In either case, the natural governing equation is a nonlinear 
advection diffusion equation that must be closed with the use of another mass balance 
equation. When considering the advection term, it is the author's opinion that equa- 



tion (2.17) is the more natural choice. The reason for this is that the bulk velocity, 

v 9 , is likely more naturally measured as compared to that of the species velocity. This 

chapter concludes our discussion on pore-scale modeling. We now turn our attention 

to building macroscale models, but in doing so we keep in mind the diffusion models 
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at the pore scale and use cues from this scale to help make proper assumptions about 
the larger scale. 



3. Hybrid Mixture Theory 

In this chapter we use a combination of classical mixture theory and rational 
thermodynamics (henceforth called Hybrid Mixture Theory (HMT)) to study novel 
extensions to Darcy's law, Fick's law, and Fourier's law in variably saturated porous 
media. This approach was pioneered by Hassanizadeh and Gray in the 70's and 80's 
[39| HU HU 05] and later extended by Bennethum, Cushman, Gray, Hassanizadeh, 
and many others [291 EH SU EI] to model multi-phase, multi-component, and multi- 
scale media. HMT involves volume averaging, or upscaling, pore-scale balance laws 
to obtain macroscale analogues. The second law of thermodynamics is then used 
to derive constitutive restrictions on these macroscale balance laws. Constitutive 
relations are particular to the medium being studied, and hence depend on a judicious 
choice of independent variables for the energy of each phase in the medium. There are 
many excellent resources for the curious reader to gain a more thorough understanding 
of HMT (eg [3TH 1ST]). For that reason we will not derive every identity along the 
way. Instead partial derivations of the identities necessary to understand the present 
application of HMT are presented. 

To begin this overview we consider the upscaling of pore-scale balance laws (con- 
servation laws) via a mixture theoretic approach. The subsequent sections in this 
and the next chapter introduce the entropy inequality and it's exploitation to derive 
constitutive laws. A judicious choice of independent variables for the energy of each 
phase in the medium is chosen and is used to derive novel versions of Darcy's, Fick's, 
and Fourier's laws. These constitutive equations will be used in subsequent chapters 
to develop models for moisture transport in variably saturated porous media. 

3.1 The Averaging Procedure 

When considering a porous medium one cannot avoid discussing the various scales 
involved. This particular work deals with two principal scales: the microscale and 
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the macroscale. At the microscale the phases are separate and distinguishable. Typ- 
ical microscale porous media will have pores that measure on the order of microns 
to millimeters (depending on the type of solid). At the macroscale the phases are 
indistinguishable and the typical measurements range from millimeters to meters. 
The macroscale is where most physical measurements are made, and as such, we seek 
to derive governing equations that hold at this scale. The microscale structure may 
vary dramatically for different media depending on the type of solid phase and the 
microscale behavior of the fluid phases. As such, the microscale geometry can have a 
dramatic influence on flow and phase interaction. 

For any given phase at the pore scale the mass, linear momentum, angular mo- 
mentum, and energy balance laws must hold. The problem is that it is difficult to 
obtain geometric information everywhere at this scale For this reason we seek to 
average (or upscale) the microscale balance laws to the macroscale. 

There are many methods for mathematically averaging balance laws. Here we 
choose the simplest method of weighted integration. Before introducing the technical 
details of the weighted integration we must first introduce the concept of a Represen- 
tative Elementary Volume and local geometry in a porous medium. This elementary 
volume will become our basic unit of volume throughout this research. The follow- 
ing discussions closely follow and paraphrase those of Bear [5], Bennethum [121 H3], 
Hassanizadeh and Gray [391 S3] , Weinstein [5Tj . and Wojciechowski [86J. 



3.1.1 The REV and Averaging 

In this work we consider unsaturated porous media. Characteristic to these media 
is the juxtaposition of liquid, solid, and gas phases within the pore matrix. We make 
the assumption that a representative elementary volume (REV), in the sense of Bear 
[5], exists at every point in space. To properly define the REV we first define the 
porosity. 
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Consider a sequence of small volumes within a porous medium, (SV)k, each with 
centroid x e M 3 . For each k, let (8V VO id)k be the volume of the void space within 
(SV)k- The porosity for the k th volume is given as the ratio 

, (8Vvoid)k 



{SV)i 



(3.1) 



Generate the sequence, {4>k}, by gradually shrinking (SV)k about x such that (SV)i > 
(SV)2 > (SV)3 > ■ ■ ■ . As k increases, the porosity will certainly fluctuate due to 
heterogeneities in the medium. As (5V)k shrinks there will be a certain value, k = k*, 
such that for k > k* the fluctuations in porosity become small and are only due to 
fluctuations in the arrangement of the solid matrix. If (5V)k is reduced well beyond 
(5V)k* the sequence of volumes will eventually converge to x. The point, x only 
lies within one phase, so the limit of the sequence of porosities will either be or 1 
(completely in the void space or completely in the solid). This indicates that there 
will be some other intermediate volume, (SV)k** < (8V)k*, where the sequence of 
porosities begins to fluctuate again as k gets larger. We define the REV, 8V, as any 
particular volume (5V)k** < SV < (5V)k*- Without loss of generality we can simply 



choose 5V = (8V)k**- Figure 3.1 illustrates two typical sequences of porosities, </>&, as 



the volume is decreased (right to left). 



Porosity 




REV 



Domain of 

homogeneity 



-> Volume 



k = k** 



k = k* 



Figure 3.1: Illustration of the definition of the REV via a sequence of porosities 
corresponding to a sequence of shrinking volumes. (Image similar to Figures 1.3.1 
and 1.3.2 in Bear 0) 



21 



Micro Scale 



REV 



Macro Scale 




Figure 3.2: Cartoon of the microscale, REV, and macroscale in a granular soil. The 
right-hand plot depicts the mixture of all phases. 



Consider now a coordinate system superimposed on the porous medium. Let x 
be the centroid of the REV, and let r be some other vector inside the REV. Define 
the vector, £, as a vector originating from the centroid of the REV such that 



r = x + £. 



(3.2) 



We can now view £ as a local coordinate in the REV as in Figure 3.3 




Figure 3.3: Local coordinates in and REV. 



Define the phase indicator function as 



la(r,t) 



1, r G a-phase 
0, r a-phase 



(3.3) 



where r is a position vector as indicated in Figure |3.3| The averaging technique 

involved multiplying a micro scale quantity (such as density) by 7 a and integrate 

over the REV. This effectively smears out the phases. A consequence of this is that 
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the averaged value may not accurately represent the actual values being measured at 
the pore scale. A further mathematical complication arises since the integrations may 
not make sense in the traditional (Riemannian) sense. Therefore, we must understand 
all of the following mathematics in the distributional sense (integrals are understood 
to be Lebesgue, and derivatives are understood to be generalized derivatives). For 
more specifics on these mathematical tools see standard graduate texts on functional 
analysis (eg. [59|). 

To find the volume of the a phase in the REV we simply integrate 7 Q over the 
REV. Define this volume as \8V a \: 



\5V a \ = / j a (r,t)dv = / -y a (x + S,t)dv(£). 

J8V JSV 



(3.4) 



The a-phase volume fraction, e a , is defined as[j 



e«(x,t)= 1 -^. (3.5) 

Since < \5V a \ < \SV\ it is clear that < e a < 1. Furthermore, since the REV is 
made up of all of the phases, 

£> = !■ (3.6) 

a 

The volume fraction is the first example of a macroscale variable. That is, it is 
a variable that describes a pore-scale property but is upscaled to the larger, more 
measurable, scale. 

It is useful to note that there are two main types of averaging that will be used: 
mass averaging and volume averaging |l3j UTJ 1ST] . Let ifrf be the j th constituent of 
some quantity of interest. To volume average ipi we define 

W = T^Tl I ^ j (r,t) la (r,t)dv(£), (3.7) 

\oV a \ Jsv 



^^Note: the notation "e a " for the volume fraction is not necessarily standard. Some authors use 
"</>"", "# Q ", or "n Q ". Furthermore, the superscript notation is sometimes replaced by subscripts. 
The present notation is chosen to be consistent with the primary references for Hybrid Mixture 
Theory mentioned in the introduction to this chapter. 
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and to mass average ipi we define 



^ = I a^iat/ I / ^{r,tW{r,t) la {r,t)dv{i). {3.1 

{p 3 ) \oV a \ Jsv 



Implicitly in (3.8) we see that density is volume averaged. That is, 

<p J 7 = Jr, / P*{r,tfta{r,t)dv(t). (3.9) 

\oV a \ Jsv 

(Note: some authors use the mass averaged notation on density even though it is 
technically volume averaged. Given the definition of a mass averaged quantity there 
is usually little confusion.) 

The basic rules of thumb for deciding whether to volume or mass average were 
originally proposed by Hassanizadeh and Grey in 1979 j^5] . They propose four cri- 
teria, listed below, for making this decision. In these criteria it is emphasized that 
the microscale quantities correspond to small scale pre-averaged quantities, while 
macroscale quantities are defined via the averaging process. 

1. "When an averaging operation involves integration, the integrand multiplied 
by the infinitesimal element of integration must be an additive quantity. For 
example, the internal energy density function, E, is not additive, but the total 
internal energy, pEdv is additive and an average defined in terms of this quantity 
will be physically meaningful." 

2. "The macroscopic quantities should exactly account for the total correspond- 
ing microscopic quantity. For example, total macroscopic momentum fluxes 
through a given boundary must be equal to the total microscopic momentum 
fluxes through that boundary." 

3. "The primitive concept of a physical quantity, as first introduced into the clas- 
sical continuum mechanics must be preserved by proper definition of the macro- 
scopic quantity. For instance, heat is a mode of transfer of energy through a 
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boundary different from work. The definition of macroscopic heat flux must 
also be a mode of energy transfer different from macroscopic work." 

4. "The averaged value of a microscopic quantity must be the same function that is 
most widely observed and measured in a field situation or in laboratory practice. 
For example, velocities measured in the field are usually mass averaged quan- 
tities; therefore, the macroscopic velocity should be a mass averaged quantity." 
This ensures applicability of the resulting equations. 

In the upscaling procedure to follow we wish to apply a weighted integration to a 
pore-scale balance law (a partial differential equation). This will involve terms such 



as 



<v 



ladv, 



isv dt 

and to ensure that we properly define the macroscale variables as either volume 

or mass averaged quantities, we need a theorem that allows for the interchange of 
integration and differentiation. This theorem is due to Whitaker and Slattery [751 E21 
183] and a generalization of this theorem is due to Cushman [29]. 

Theorem 3.1 (Averaging Theorem) Ifw a p is the microscopic velocity of the aft 
interface and n a is the outward unit normal vector of5V a indicating that the integrand 
should be evaluated in the limit as the a(3 interface is approached from the a side, 
then 

i r df 



\sv\ 



sv 



dt 



d_ 
dt 



J a dv(£) 



1 



\SV\ 



5V 



\_W\J8V 

Vf la dv(£) 

1 



fl*dv(£) 



^ \ by \ J** 



(3.10a) 



fladv(£) 
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where f is the quantity to be averaged. 
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Keep in mind that / could be a scalar or a vector quantity. In the latter case, the 
symbol / is replaced with / and appropriate tensor contractions are inserted. 
The averaging procedure is now carried out in the following steps: 

1. State the pore-scale balance law for a particular species (or phase). 

2. Multiply the equation by j a . 

3. Average each term over the REV. 



4. Apply Theorem 3.1 to arrive at terms representing macroscale quantities. 



5. Define physically meaningful macroscopic quantities. 

We now turn our attention to averaging pore-scale balance laws in the sense listed 
above. In the following discussion, v^ is the microscopic velocity of constituent j, w a p. 
is the velocity of the j th constituent in the a/3 interface, and n a is the outward unit 



normal vector of SV a . A full nomenclature index can be found in Appendix B.l 



3.2 Macroscale Balance Laws 

As it is the simplest balance law, let us first consider the mass balance equation 

for a single constituent: 

Dp> ... 

-^- + p 3 V-V 1 =p 3 r 3 . (3.11) 

Here, p- 7 is the density of the constituent, v^ is the velocity of the constituent, and 
any source of mass from chemical reactions between the constituents is given as r J . 
Recall that the material (Lagrangian) derivative is 

*£-%+*■* <•> 

This derivative contains the usual Eulerian derivative along with an advective term. 



Written in terms of the Eulerian time derivative, (3.11) is 

dp? 



+ V- (pV) = pPf i . (3.13) 
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While this is specifically the mass balance equation, it takes the prototypical form of 
all balance laws: a time derivative plus a flux is equal to any source. 

The constituent momentum balance can be written in a similar manner: 



Ot 



V ■ [p'v-' ® v^ — t J ') = (P (g + i + r^v 



(3.14) 



Here, V is the Cauchy stress tensor on species j, and the sources on the right-hand 
side are gravity, momentum transfer from other constituents, and momentum gained 
from chemical reactions respectively. These equations describe the change in mass 
and momentum over time and space within a specific constituent. They are sufficient 
field equations for modeling systems composed of a single phase gas, liquid, or solid, 



but equations (3.13) and (3.14) are insufficient for modeling multiphase and multi- 



constituent systems as a mixture because the interactions between the phases and 
constituents are not present. 

In this work we consider a porous medium consisting of a solid phase and two 
fluid phases with multiple constituents within each phase. The phases will be denoted 
as a = l,g, and s for liquid, gas, and solid respectively. The constituents will be 
enumerated j — 1 : N (using MATLAB-style notation to indicate j = 1, 2, ... , N). 
The following derivations follow similar derivations given by Gray |43j . Weinstein [81] 
and Wojciechowski [86J. 

3.2.1 Macroscale Mass Balance 

To obtain macroscale equations in multiphase and multi-constituent media we 
multiply a constituent balance equations by the phase indicator function, j a , integrate 



over 8V, and divide by \SV\. Applying the averaging theorem (3.1 ) to the appropriate 



terms in equation (|3.13|) we have 
1 



\8V\ 



5\- 
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d_ 

d~t 



(s a P A - V -*- / fPwcf, ■ n a da } (3.15a) 
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\SV\ 
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fpf^dv = s a pi fi . 



(3.15b) 
(3.15c) 



a\ JSV 



Substituting equations (3.15)(a)-(c) into equation (3.13), recognizing the volume 
fraction terms, and recognizing the averaged mass and velocity gives the upscaled 
mass balance equation: 

d (e a p~i C 



Ot 



V • £ a pi vi 



^ a \W\J5A aP 



(3.16) 



Rewriting equation (3.16) in terms of the material time derivative and defining 






-Q: 



e« J = 



\6V C 



a\ JSA a p 

■a 



P 3 ( w a/3j — v J ) ■ n a da, and 



(3.17) 
(3.18) 
(3.19) 

(3.20) 



respectively to be the averaged mass over 5V a , the mass averaged velocity, the net 
rate of mass gained by constituent j in phase a from phase /3, and the rate of mass 
gain due to interaction with other species within phase a, we get 



^^W'v-^E^ 



f aj . 



(3.21) 
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Next we define the corresponding bulk phase variables so that the macroscale 
equations are consistent with experimentally measured terms as much as possible. 
Define: 



TV 



p a = j~y J; and 



c°y 



7" 



(3.22) 



(3.23) 



respectively to be the mass density of the a phase, and the mass concentration of the 
jth cons tituent in the a phase. If equation (3.21) is rewritten as 

1 p ' + v • ( £ y^^) = J2 e? + ^ 



fca 



and then summed over j = 1 : N we obtain a mass balance equation for the a phase: 

d (e a p c 



dt 



+ V ■ (e a p a v c 



E 



e /3- 



(3.24) 



^JV ~ctj 



Here we have used the fact that eg = Yuj=t &b '■> ^ e ra ^ e °f mass transfer to the a 
phase from the j3 phase is the sum of the rates of mass transfer to each individual 
constituent in the a phase from the /3 phase. 

Now use the definition of the material time derivative to write the mass balance 
equation for the a phase as 

D Q (e a p c 



Dt 



+ e p V -v = 2^ e /3' 



(3.25) 



where the following restrictions have been applied: 

JV 

J2r aj =0, Va, and 



i=i 



a j3^a 



0, j = l:N. 



(3.26) 
(3.27) 



Restriction (3.26) states that the rate of net gain of mass within species a from 



chemical reactions alone must be zero. Equation (3.27) states that the rate of mass 

gained by phase a from phase /3 is equal to the rate of mass gained by phase /3 from 

phase a. 
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3.2.2 Macroscale Momentum Balance 



We now turn our attention to the momentum balance equation, (3.14). We can 
apply the same principles to upscale this equation (for full details see [SI])- The 
macroscopic linear momentum balance equation for constituent j in the a phase is 

£ V j ^^ - V ■ (e a f J ) - e a p a ig a i = T 3 + ^ T^ , (3.28) 

and the macroscopic linear momentum balance equation for the a phase is 

EapaD ~& ~ V ' ^ " ^^ = 5-^' (3 ' 29) 

where £"■?' and t a are the Cauchy stress tensors for the species and the phase and 
To and T a are momentum transfer terms. Most specifically, for the momentum 
transfer terms, the former represents momentum transfered to constituent j in the 
a phase through mechanical interactions from phase (3, and the latter represents 
the momentum transfered to phase a through mechanical interactions from phase j3. 
Also notable is the i 3 term. This term represents the rate of momentum gain due 
to mechanical interactions with other species within the same phase. 

In the processes of deriving these equations the following restrictions were en- 
forced: 

N 

Y^ (P +r a 'v aj ' a ) = Va, and (3.30) 

EE(^ + ^?) =0 j = l:N. (3.31) 



Restriction (3.30) states that linear momentum can only be lost due to interactions 



with other phases (not within the species), and restriction (3.31) states that the 



interface can hold no linear momentum. The comma in the superscript of (3.30) 



indicates a relative term: v a i' a = v aj — v a . For a complete list of notation see 



Appendix B.l 
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Lastly, to tie the momentum transfer and stress tensor for the a phase to those 
of the species we note two identities that were used in the derivation: 

N 

I" = Yl (f ' ~ P a3 v a " a ® v a "' a ) (3.32) 

3=1 

N 

f; = J2{T a /+e^v^). (3.33) 

i=i 

These identities will be used later and so are presented here for conciseness. 

3.2.3 Macroscale Energy Balance 

The derivations for the macroscale angular momentum and energy balance laws 
are more algebraically complicated. The angular momentum equation will not be 
used in this work since we assume that we're dealing with granular-type media where 
the angular momentum balance results in the solid phase Cauchy stress tensor being 
symmetric [151 H31 EI]- The energy balance equation, on the other hand, will allow 
us to derive a novel form of the heat equation in porous media. For this reason we 
state the full equation here. 

Applying the same routine as in the mass and linear momentum equations we 
arrive (after significant simplification) at a balance law for the energy in species j: 

£ a paj D a >{e a >) _ v _ ( £ a q a^ _ £ ap . Vv * 3 _ £ « p « 3 ^ = ga 3 + £oy ^ g^ 

(see [3 E] for details on the derivation). Here, h aj is the external supply of energy, 
e aj is the energy density, q aj is the partial heat flux vector for the j th component 
of the a phase, Q aj is the rate of energy gain due to interaction with other species 
within the a phase, and QJ is the rate of energy transfer from the phase to the a 
phase not due to mass or momentum transfer. 

Again, following the derivation of [81], the bulk phase energy equation is 

e a p aL ^f- - V • (e a q a ) - e a f : Vv a - e a p a h a = J^ Q% ( 3 - 35 ) 

J—J L 
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where 

To arrive at this form of the energy equation we enforced the following restrictions: 



N 

E 

EE 



()"' + C 3 v a ^ a + r a > ( e a > + - (v a >' a f 



Q\ 



T?-v a > + < 



3 I p"3 _|_ 



(v a J' a ) 2 



Va, and 



j = 1: N. 



(3.36a) 
(3.36b) 



Restriction (3.36a) states that energy gained or lost due to species interactions within 



the a phase must be gained or lost due to interactions with other phases. Restriction 



(3.36b) states that the rate of energy gained or lost by one component in one phase 



must go to another component or phase. That is, this second restriction states that 
the interface retains no energy. 

A system of equations governed by mass, momentum, and energy balance requires 
each of the upscaled equations listed. A count of the variables indicates that there 
are far more variables than equations. It is at this point where we need a method 
for deriving constitutive equations for these remaining variables. The method cho- 
sen for this work uses another macroscale balance law based on the second law of 
thermodynamics. 

3.3 The Entropy Inequality 

The development of constitutive laws is central to the modeling process. As 
we mentioned previously, this has historically been a process of fitting mathemati- 
cal models to empirical evidence. The construct of Hybrid Mixture Theory (HMT) 
couples the averaging theorems discussed in the previous section and the second law 
of thermodynamics to provide us with restriction on the form of the constitutive 
relations; hence narrowing down the experiments required to those that are thermo- 
dynamically admissible. It is then up to the experimentalists to verify and refine these 
models. Both theoretical and experimental directions of study have their merits, but 
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putting the constitutive equations on a firm theoretical footing is ultimately preferred 
whether it is before or after the experiments are run. In this section we give a brief 
derivation of the upscaled entropy inequality, and we then use this inequality, along 
with a judicious choice of variables, to derive constitutive equations for unsaturated 
porous media. 

3.3.1 A Brief Derivation of the Entropy Inequality 

The second law of thermodynamics states that entropy will never decrease as a 

system evolves toward equilibrium [U [2T] . The microscale entropy balance equation 
that describes this phenomenon is 

p> ^f- + V • 4? - p>V = ff + fV + A, (3.37) 

where rf is the entropy density of constituent j, (p 3 is the entropy flux, V is the 
external supply of entropy, ff is entropy gained from other constituents, and A is the 
entropy production. Since the second law of thermodynamics must hold we know 
that A > for all time. 



Applying Theorem 3.1 to equation (3.37) and defining appropriate macroscale 



definitions of the variables gives the upscaled entropy balance equation: 

e a p a > ^pL _ V ■ {e a 4> a >) - e a p a ib a > = ^ $^ + fj aj + k aj , (3.38) 



where the terms on the right-hand side of equation (3.38 ) represent transfer of entropy 
through mechanical interaction, entropy gained due to interactions with other species, 
and the rate of entropy generation respectively. 

Next, assume that the material we are modeling is simple in the sense of Coleman 
and Noll [27]. This means that we assume that the entropy flux and external supply 



are due to heat fluxes and sources respectively. To remove the dependence on external 



heat sources we add (1/T times) the upscaled conservation of energy equation (3.34) 



\a, 



e a p a > D n 6 3 - V ■ (e a q a i) - e a fi : Vv aj - e a p^h a ^ = Q a > + ^ Q c p . 
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At this point we perform a Legendre transformation in order to convert convert 
internal energy, e a \ to Helmholtz potential, ip aj (see any thermodynamics text, eg 

my- 

^i = e *i - jfiT. (3.39) 

This is done because internal energy has entropy as a natural independent vari- 
able, and entropy is difficult to measure experimentally. It should be noted that 
the Helmholtz potential is only one choice of thermodynamic potential we could have 
made. This is done primarily for historical reasons, but the Gibbs potential and possi- 
bly the Grand Canonical potential could have also been viable choices. The appeal of 
the Helmholtz potential is that it naturally has independent variables of temperature 
and volume (or, in intensive variables, density). 

To arrive at a simplified entropy inequality for the total production of entropy 
(across all constituents and phases) we now solve for A^ and then sum over a — l,g,s 
and j — 1 : N. This step requires significant algebra so the details of the derivation are 
omitted for brevity sake. After much simplification, the entropy inequality becomes 
e a p a (D a ^) a a D s T 



Z^\ T \ Dt ' Dt 
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r'\[^ a +^V 



a,s „,a,s 



(3.40) 



where A is the rate of entropy generation. 



Several new terms have appeared in (3.40). First, d a = CVv a ) sym is the rate of 



deformation tensor (also known as the strain rate). As before, terms with a comma 
in the superscript are relative terms: v a,b = v a — v b . 



Several identities were needed to derive (3.40). A complete list of these identities 



has been included in Appendix |B.3| The next step is to expand the Helmholtz 
potential in terms of constitutive independent variables that describe our system. 
This allows freedom to make choices about which variables control behavior of the 
system. The choice of these variables is generally non-trivial so in the next section 
we discuss motivations for the choice of variables. 
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4. New Independent Variables and Exploitation of the Entropy 
Inequality 

Now that we have an expression for the entropy inequality we must choose a set 
of independent variables that describes our system of interest. We seek to describe 
a multiphase system where the solid phase may undergo finite deformation, where 
the relative saturations of the two fluid phases vary in time and space, and where 
phase changes between the fluids possibly occurs throughout the porous medium. 
Hassanizadeh and Gray have modeled similar media in the past [4"T1 HU H5j. These 
models include effects from common interfaces, common lines (where three phases 
meet), and common points (where four phases meet). These models are very thorough 
and follow the same HMT approach. The down sides to their models, in the author's 
opinion, are three fold: (1) the complexity of the resulting equations is such that in 
order to use these equations a host of simplifying assumptions must be made, (2) the 
thermodynamics of the common points and lines make sense physically but are likely 
negligible relative to other effects, and (3) constitutive equations must be derived for 
transfer rates between interfaces, common lines, common points, and phases. This 
final drawback indicates that a detailed knowledge of the pore-scale physics must 
be somehow upscaled. Approaches have been taken recently to do just this, but the 
proposed theories have not yet gained widespread acceptance. Examples of such work 
include those of Gray et al. [I2J HO] 

In the present approach we choose not to directly model interfaces and instead 
strive to eventually write our governing equations in terms of the macroscale chemical 
potential. The chemical potential is known from physical chemistry and thermody- 
namics as a generalized driving force that is a function of pressure and temperature. 
It is well known that mass transfer from liquid to gas states is driven by gradients 
of chemical potential [SB], so if we can write constitutive equations (such as Darcy's, 
Fick's, and Fourier's laws) in terms of this potential we can possibly couple the rele- 
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vant effects into much simpler governing equations; for example, equations that track 
changes in chemical potential instead of pressure or concentration. The immediate 
drawback to the present modeling efforts is that the recent work by Hassanizadeh et 
al. seems to indicate that saturation and capillary pressure are linked to the amount 
of interfacial area between phases within the medium [20j EH] • In the present work 
we will not directly model the fluid-fluid and fluid-solid interfaces. We proceed with 
the present modeling effort despite the results proposed by Hassanizadeh et al. We 
will discuss this drawback as we run up against it in future sections and chapters. 

4.1 A Choice of Independent Variables 

In this section we present a choice of independent variables for the Helmholtz free 
energy (potential) so as to expand the entropy inequality and to derive the relevant 
forms of Darcy's law, Fick's law, and Fourier's law. These variables are known as 
constitutive independent variables as they represent a postulation of the variables that 
control the energy in the system. "Deriving physically meaningful results depends on 
our ability to relate thermodynamically defined variables to physically interpretable 
quantities" |81j . To that end, we use our a priori knowledge of thermodynamics to 
choose some of the variables. For the remainder of this work we restrict our attention 
to a three-phase system consisting of an elastic solid, a viscous liquid phase, and a gas 
phase. To begin the modeling process we assume that each of these phases consists 
of N constituents (also called species or components), and all interfacial effects are 
neglected. Examples of the constituents include dissolved minerals in the liquid, 
species evaporated into the gas, or precipitated minerals associated with the solid 
phase. 

The motivation for choosing some of the variables is relatively trivial. For ex- 
ample, to allow for a heat conducting medium, temperature, T, and the gradient of 
temperature, VT, are included in the list of independent variables. The pore space is 

expected to be variably saturated with the two fluid phases so the volume fractions, 
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e l and e 9 , must be included in the set of variables. The fact that J2 a £ Q = 1 precludes 
us from using all three volume fractions since they are not independent of each other. 
In future chapters we will further restrict this assumption since for a rigid solid phase 
the sum of the fluid phase volume fractions is equal to the fixed porosity 

e l + e 9 = e. (4.1) 

The reason for not making this assumption initially is that it allows us to develop 
models for deformable media as well as for media with a rigid solid phase (hence, a 
more general model may be derived from these assumptions later if necessary). 

Recall from thermodynamics that the change in extensive Helmholtz potential, 
A, with respect to volume is minus the pressure: dA/dV = —p. In terms of intensive 
variables this means that p 2 dip/dp = —p. To remain consistent with the extensive 
definition of the Helmholtz potential, the densities must then be included in the set 
of independent variables. Given the fact that there are N constituents in each phase, 
this could be done in two different ways: (1) we could include the mass concentrations, 
C"- 7 , for j — 1 : N — 1 along with the phase density, or (2) we could include all of 
the constituent densities, p aj for j — 1 : N, Bennethum, Murad, and Cushman 
[13] . and also Weinstein [81] took the first of these options when using HMT to 
derive constitutive relations involving chemical potentials. The trouble with this 
approach is that the mass concentration of the N th constituent is dependent on the 
mass concentrations of the previous N — 1 constituents (since the concentrations 
sum to 1). These results indicate that the behavior of the constituents depends 
on how they are labeled instead of simply being independent. Various techniques 
were successfully developed in [15] to deal with this complication. To avoid these 
complications we choose the second option and include the species densities, p aj for 
j — 1 : N. Since each constituent is free to move within each phase, the spatial 
gradients of the species densities, Vp' J and Vp* , are also included. 
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Darcy's law and Fick's law are classical empirical expressions for creeping flow 
and constitutive diffusion. Darcy's law is a statement about the relative velocity of 
a fluid phase in a porous medium, and Fick's law is a statement about the relative 
diffusive velocity of a species within a phase. Since we seek novel forms of these two 
laws we include v a,s and v a i' a for a = l,g, s in the list of independent variables. It 
should be noted that neither of these variables is objective in the sense that they are 
not frame invariant. This poses a problem since any governing equation should not 
depend on an observer's frame of reference. In [31] , Eringen proposed a modification 
to Darcy's law that creates a frame invariant relative velocity. The new terms needed 
for this new relative velocity are second order and are assumed to be negligible in 
Darcy flow. A similar argument can be used for Fick's law. 

The reasoning given in the previous few paragraphs leads us to the set of inde- 
pendent variables for ip a to include: 

T, V2V,£ 9 ,p' J ,P 9j , VA V(P,v l > 8 ,v 9 > 8 ,v li '\v 9i ' 9 , and v Sj ' s 

where j — 1 : N. It is apparent, now, that solid-phase terms corresponding to the 
density and gradient of density are missing. The principle of equipresence, from 
constitutive theory in continuum mechanics, states that "all constitutive variables 
are a function of the same set of independent variables" jHH]- To give symmetry 
between the phases we include p s and Vp s . The Stokes assumption for the Cauchy 
stress tensor in a viscous fluid states that stress is the sum of the fluid pressure and 
the strain rate. For this reason we include the strain rate (also known as the rate 
of deformation tensor) for the fluid phases: d and d 9 . The theory of equipresence 
also states that if we include strain rate in the fluid phases then we must include a 
comparable term in the solid phase. 

A natural choice of variables for the solid phase are the solid phase volume frac- 
tion, density, and the (averaged) strain. Weinstein [Hi] pointed out that these three 

variables are not independent, as explained below, and used a modified set of inde- 
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pendent variables for the solid phase. The same modified set will be used here, so 
the following simply states Weinstein's results with brief derivations. 

Let J s be the Jacobian of the solid phase given by J s =det((-F s ) T • -F s ), where 
-F s is the deformation gradient 

vs _ - /, ( 42 ) 



= dX^ 

x s is the Eulerian coordinate, and X s is the Lagrangian coordinate. Using standard 
identities from Continuum Mechanics, the Jacobian can be rewritten as 

J s = det (22| s + I) . (4.3) 

Furthermore, through the conservation of mass, the Jacobian is also a scaling factor 
for volumetric changes, J s = (SqpI)/(s s p s ). This clearly shows the dependence of the 
three variables. To mitigate this issue, Weinstein [HI] adopted ideas from solid me- 
chanics and considered a "multiplicative decomposition" of the deformation gradient, 
F_ s , and the Green's deformation tensor, C_ s , as 

Q s = (J s ) 2/3 a s , (4.4) 

where (J S ) 1//3 J and (J s ) 2 / 3 / represent volumetric deformation, and U F_ and G are 
the modified deformation gradient and the modified right Cauchy-Green tensor, re- 
spectively." With this modification to the solid strain, the solid phase variables we 
consider here are J s ,d , C Sk and VC Sfc where k — 1 : N — 1. We note here that 
in order to get physically meaningful results for phase change, we include the same 
components, so that C Sj ,C l: >, and C 9j all refer to the same component. Pairing the 
mass concentrations and the Jacobian gives a description of the density of the solid 
phase, and the modified Cauchy-Green tensor is used in place of strain. 

The principle of equipresence states that all of the constitutive variables must be 

a function of the same set of the postulated independent variables. In particular, we 
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postulate that the Helmholtz potential for each phase is a function of the following 
set of variables: 

{T, VT, e\ e\ p l \ff\ Vp l >, Vp 9 ', v a >%£, § , v a >' a , J S ,Q S , C s \VC Sk } , (4.6) 

where a = I, g, s; j — 1 : N and k — 1 : N— 1. We postulate that a three phase porous 
medium with an elastic solid phase and N constituents per phase can be modeled by 



set (4.6). 



4.1.1 The Expanded Entropy Inequality 



Consider now that the first line of the entropy inequality, (3.40), contains a ma- 
terial time derivative of the Helmholtz potential for the a phase. Using the identity 

D a (-) D s (-) as _.. 

+ v a ' s -V(-), (4.7) 



Dt Dt 

and applying the chain rule, the entropy inequality can be expanded to include each 
of our constitutive independent variables. The central idea to the exploitation of the 
second law of thermodynamics is that no term in the entropy inequality can take 
values such that entropy generation is negative. A close examination of the expanded 
entropy inequality reveals that there are many terms that show up linearly. In these 
linear terms we notice some that are neither independent nor constitutive. Examples 
of such coefficients are VT, VC% V//\ Vp*',^,#, v l ' s ,v 9 ' s ,v s ^' s , v l *'\ v 9 ^ 9 , T, p a \ 
and Vi> aj ' a (where the dot notation (e.g. T) indicates a material time derivative). 
Loosely speaking, we have no control over these variables and they could take values 
that violate the second law. For example, take as a thought experiment a process 
where all of these variables except T are zero. From Bennethum [ID] , 



"Since none of the other terms in the entropy inequality are a function 
of T, by varying the value of T we can make the left-hand side of the 
entropy inequality as large positive or as large negative as we want - hence 
violating the entropy inequality. Since the entropy inequality must hold 
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for all processes (including those for which T is any value), the entropy 
inequality can be violated unless the coefficient of T is zero." 



In order not to violate the inequality in (4.13), the coefficients of all of these factors 
must be zero. This implies that terms such as Y^ Q ( £a p a ^r) are zero anc ^ w ^ 



therefore be left out of the expansion of (3.40) for brevity. The time rates of change 



of volume fractions are not this type of variable since they are constitutive; that is, 
we assume a rule for the time rates of change of volume fractions that depends on 
the specific medium of interest. 



With this simplification in mind, (3.40) becomes 
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The next step is to enforce two additional relationships using Lagrange multipli- 
ers. In doing so, the Lagrange multipliers become unknowns of the system. We will 
see in subsequent sections that the Lagrange multipliers are associated with partial 
pressures and chemical potentials of species in the fluid phases. The first relationship 
considered is the dependence of the diffusive velocities: 

N 

J2{C aj v aj ' a ) = 0. (4.9) 

i=i 

One can see this since 

N N N 

J2 C a >v a >' a = J2 C aj v aj - J2 ° aJva = v a -v a = 0. (4.10) 

j=i j=i j=i 

The implication is that if we know the concentrations and diffusive velocities of the 

first N — 1 constituents, then we would know the concentration and diffusive velocity 

of the N th constituent. Multiplying by the density, taking the gradient, and using 

the product rule gives the following relationship: 

(N \ N 

Y^fptv 01 ** I = ^ {p aj Vv a *' a + v^Vff 1 *) = 0. (4.11) 

3=1 / 3=1 

Following Bennethum, Murad, and Cushman [15J, we enforce this relationship with 
a Lagrange multiplier so as to account for the N th term dependence. 
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The second relationship to be enforced with Lagrange multipliers is the mass 



balance equation for each of the constituents (3.21): 
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+ e a p^V ■ v a > = E e? + r a >- 



Let A o[d denote A from equation (4.8), and let X aj and X aN be the Lagrange 



multipliers for the mass balance and N th term dependencies, (4.11), respectively. 
The entropy inequality is rewritten as follows: 
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After a significant amount of algebraic simplification (with no additional physical 
assumptions), this yields the following form of the entropy inequality: 
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The exploitation of equation (4.13) will be the source of all of the constitutive 



relations for the remainder of this work. The next section outlines the details of this 

exploitation to form constitutive relations specific to multiphase media governed by 
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our choice of constitutive independent variables, (4.6). 



4.2 Exploiting the Entropy Inequality 



In this section we exploit the entropy inequality, (4.13), in the sense of Colman 
and Noll [27] . The basic principle here is that, according to the second law of ther- 
modynamics, entropy is always non-decreasing as time evolves. This fact is used to 
extract constitutive relationships from the entropy inequality. Not every result from 
this exploitation is relevant to the current study, so we only present the more notable 



and useful results in the next subsections. Furthermore, we exploit equation (4.13) 
with an eye toward deformable, multiphase, media. The assumption of deformable 
media will be removed in the future, but this leaves open the possibility of returning 
to these results for future work. For an abstract summary of how the exploitation 
of the entropy inequality works, along with subtle but important assumptions, see 
Appendix [Cj 

4.2.1 Results That Hold For All Time 



As mentioned in Section 4.1.1, several of the terms that appear linearly in the 
entropy inequality have factors that are neither independent nor constitutive. We now 
use this fact to derive relationships that must hold for all time in order to not violate 
the second law of thermodynamics. To illustrate this point consider the coefficient 
of T. If this coefficient is set to zero we recover with the thermodynamic constraint 
that temperature and entropy are conjugate variables, 

dib a 

-^r = ~V a . (4-14) 

This is a classical result known from thermodynamics. 
4.2.1.1 Fluid Lagrange Multipliers 
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For the gas and liquid phases, the definitions of the Lagrange multipliers stem 
from the coefficient of p aj and Vi> Qj ' Q . Setting the coefficient of p aj to zero gives the 
definition of the Lagrange multiplier for the mass balance equations: 

a 

Setting the coefficient of Vt) ai,a to zero, summing over j = 1 : N, and solving for 
A Qjv yields an expression for the other Lagrange multiplier: 

1 - 

A°^ = — - Yl ii a ' + (p aj A ° J ) £ + ^ a L ( 4 - 16 ) 

3=1 

4.2.1.2 Solid Phase Identities 

Several identities for the solid phase can be derived from the terms associated with 
the time derivatives of the solid phase Jacobian, J s , and the modified Cauchy-Green, 

• s 

C terms. From the J s term we see that 

N 

3 



i=i 
Next, consider the identity 

N 



o' = 1 rv \ / 



r = 5Z (r j + p ajv ° j ' a ® ^' a ) ( 4 - 18 ) 

i=i 



resulting from upscaling the momentum balance equation. Taking the trace of (4.18) 



neglecting the diffusive terms, and substituting this into (4.17) gives a definition for 
the solid phase pressure: 

This is a generalization of the solid phase pressure found by Weinstein for saturated 
porous media in [ST] . 

The coefficient of the G term gives a relationship for the stress in the solid phase. 

This will give a generalization of the solid phase stress [SJ [9] and closely follows the 
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derivations of Bennethum [9] and Weinstein [81]. Setting the coefficient of the Q 
term to zero, left multiplying by the modified deformation gradient, F_ , and right 
multiplying by the transpose of the deformation gradient gives a relationship that 
defines the Lagrange multiplier for the solid phase, X Sj : 
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Using identity (4.18) in the stress term of (4.20), neglecting the diffusive velocities 



taking one-third the trace of the result, and using equation (4.19) for the solid-phase 



pressure yields a relationship for the solid phase Lagrange multiplier: 
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Substituting (4.21) back into (4.20) gives the following relation for the solid phase 

stress: 
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This can be rewritten as 
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(4.24a) 
(4.24b) 
(4.24c) 



3 r dQ 
The stresses above are termed the effective stress, hydrating stress for the liquid 



phase, and hydrating stress for the gas phase respectively. Equation (4.22) states 
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that the stress in the solid phase can be decomposed into the solid pressure and 
stresses felt due to the presence of the fluid phases. It is here that the modifications 
of the deformation gradient and Cauchy-Green tensors become clear. If we take the 
trace of the stress tensor then we see that (l/3)tr(t s ) = —p s ; which is how the 
solid phase pressure is measured. Therefore, this thermodynamic definition of p s is 
consistent with experimental measure. Furthermore, the effective stress and hydrating 
stresses are terms associated with the interaction between the solid and the fluids. 
For saturated porous media, Bennethum jH] states the following: 

"The effective stress tensor is the stress of the solid phase due to the 
strain of the porous matrix, and the hydrating stress tensor is the stress 
the liquid phase supports due to the strain of the solid matrix (which 
would be negligible if the liquid and solid phase were not interactive, but 
which becomes significant for swelling porous materials)." 

One final note on the solid phase stress is that the total stress in the porous 
medium is related to the pressures in all three phases. This can be seen by taking 
the weighted sum of the stresses in each phase: 

i = Yl e T = - £S p s l + e % + £l (i h + i) + e9 (§ + 1 9 ) • ( 4 - 25 ) 

a 

Taking one-third the trace of the total stress, and recalling that the effective and 
hydrating stresses are trace free, gives 

Q) tr (i) = - £s p s + £ ^ tr & + j tr (i 9 ) ■ ( 4 - 26 ) 

In order to fully understand the stresses in the fluid phases we must continue our 



examination of the results coming from the entropy inequality. Equation (4.26) is 
similar to the Terzaghi stress principle; suggesting that the fluid phases help to sup- 
port the pore space in the medium. 

4.2.2 Equilibrium Results 
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There are several more relationships that we can extract from the entropy in- 
equality. In particular, we now seek relationships between the Lagrange multipliers 
and the fluid-phase pressures. We also seek relationships for the momentum and en- 
ergy exchange terms. At equilibrium the production of entropy is minimized. Since 
this is a minimum, the gradient of A with respect to the set of independent variables 



(4.6) is zero. This indicates that the coefficients of the independent variables that 
appear linearly in the entropy inequality are zero at equilibrium. In the case of a three 
phase porous medium of this nature, we define equilibrium to be when a subset of 
the independent variables are zero. In particular, equilibrium is defined when no heat 
conduction occurs, VT = 0, the strain rates in the fluid phases are zero, dr = 0, and 
all relative velocities are zero, v a :>' a = and v a,s = 0. This definition of equilibrium 
is particular to this type of media and is chosen as it gives physically relevant and 
meaningful results. Another way to look at this is to say that equilibrium is exactly 
the state when all of these variables are zero. 



4.2.2.1 Fluid Stress Tensor 

The first notable equilibrium result comes from the coefficient of the rate of 
deformation tensor, dr. Setting the coefficient of dr to zero, eliminating the sum of 
the constituent stress tensors using the identity 

N 
3=1 

and noting that at equilibrium the diffusive velocity is zero, yields 

N 



J2 &{?> = -~tr (f) = pf>. (4.28) 



3=1 



Equation (4.28) links the Lagrange multipliers to the equilibrium pressure of the fluid 



phases. This is the classical definition of pressure in a fluid: minus one-third the trace 



of the stress tensor. Using equation (4.15) the pressure in the fluid phases can now 
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be written as 



j=l a 



With the definition of pressure in equation (4.29) we note that the coefficient of 
e 13 can now be rewritten as 

The second term is the change in energy with respect to volumetric changes, and 
is therefore interpreted as the relative affinity for one phase to another. That is, 
this term is related to the wetability of the a and solid phases by the /3 fluid phase. 
The time rate of change of volume fluid phase volume fraction is an equation of 



state (that is not yet known), but rewriting the coefficient as in (4.30) hints at the 
fact that the equation of state is related to the pressure and the wettability of the 
phases. Furthermore, pressure, wettability, and surface tension are related to capillary 
pressure; hence indicating that the equation of state for the time rate of change of 
volume fraction is related to capillary pressure. It is here that we note the drawback 
to the present modeling effort. Recall that in the present expansion of the entropy 
inequality we do not include interfacial effects. If we were to include these effects then 
a surface tension term would appear here (as shown in Hassanizadeh and Gray |41j ) 
and these terms together would more readily be associated with capillary pressure. 
More discussion will be dedicated to the exact equation of state for the time rate of 



change of volume fraction after a discussion on cross coupling pressures in Section 4^3 
and capillary pressure in Chapters [5] and [7j 

4.2.2.2 Momentum Transfer Between Phases 



The next notable equilibrium result we can extract from (4.13) comes from the 



coefficient of the fluid phase relative velocities, v^ ,s . Setting this coefficient to zero, 

recalling that VT = at equilibrium, using the definition of the fluid phase pressure, 
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(4.29), the definition of the fluid phase Lagrange multipliers, (4.15), and solving for 



the momentum transfer terms gives 
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where 7 is the other fluid phase not equal to 0. This particular result will be coupled 
with the conservation of momentum to yield novel forms of Darcy's law in Section 

axu 

4.2.2.3 Momentum Transfer Between Species 

Another notable equilibrium results comes from the coefficient of the diffusive 



velocity, v aj ' a . Equation (4.28) indicates that at equilibrium the definition of the 



Lagrange multiplier, equation (4.16), simplifies to 



1 aw _ „/," 



4> a L 



(4.32) 



This implies that, at equilibrium, the stress tensor for constituent j (from the coeffi- 
cient of Vv aj,a ) can be written as 



t a i = -E a p a i {\ a > - iff*) I - E a p a iij) a I_. 



Consider the diffusive velocity term in the entropy inequality: 
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- A Qj V (e a p a i) - A ajv : V (e>°0 \ ' v aj ' a . 

At equilibrium the diffusive velocity is assumed to be zero. By the logic used herein 
for the exploiting the entropy inequality, and given the fact that A QJV = ip a I at 
equilibrium, we observe that for each j, 

J2 T% + T* = - V (e a p a ^ a J) + \ a >V (eVO +fV (e a p a ^) . (4.34) 

This is an expression for the momentum transfer for species j in the a phase. This 
result will be coupled with the constituent conservation of momentum equation to 
derive a form of Fick's law in Section 14.5.31 

4.2.2.4 Partial Heat Flux 

To conclude the equilibrium results we examine the VT term in the entropy 
inequality. We have assumed that VT = and v a i' a = at equilibrium, so by the 
logic used above we see that the coefficient of VT must be zero and 

J2 ea q a = ° ( 4 - 35 ) 

a 

at equilibrium. This is the partial heat flux of the entire porous media and will be 



used in Section 4.5.4 to derive a generalized Fourier's law. 



4.2.3 Near Equilibrium Results 

The next step in exploiting the entropy inequality is to derive near equilibrium 
results. These results arise by linearizing the equilibrium results about the equilibrium 
state. The linearization process is simply the first-order terms of the Taylor series, 
but one must keep in mind that each of the derivatives is a function of all of the 
constitutive independent variables that are not zero at equilibrium. For example, 
if / — feq at equilibrium, then near equilibrium, / w f near = f eq + (df/d(VT)) ■ 
VT + • • • + (df /dv l,s ) ■ v l,s . This full expansion may yield terms that are not readily 
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physically interpretable. For this reason, considerable efforts must be made to relate 
the linearization constants to measurable parameters. For a thorough explanation of 
the linearization process with the entropy inequality see Appendix [Cj 



For the momentum transfer in the fluid phases, the linearization of equation (4.31 ) 
can be simply written as 

£#) =fe^) -{effg-vf*. (4.36) 

\a£fl J near W/3 / e g 

The linearization constant, Rr , is related to the resistivity of a porous medium; the 
inverse of the hydraulic conductivity. It should be noted that we have only expanded 
about one of the possible variables: v^ ,s . Strictly speaking this is incorrect and 
we should expand about all other variables which are zero at equilibrium. A more 
thorough expansion is 



E^l = £« -^rg-^ s 

W/3 / eq 

+ H? -VT + J 13 ■ v 0j,fi + LP : dP + • • • , (4.37) 



^P / n.eq W£ , eq 



where Hf and Js are second-order tensors and If is a third-order tensor. The ellipses 
at the end of this equation indicates that there are higher order terms that are not 



being written explicitly. The left-hand side of (4.37) is the rate of momentum transfer 
due to mechanical means. It is reasonable to think that this transfer term might be a 
function of fluid velocity, but the effects due to thermal gradients, diffusive velocity, 
and velocity gradients are likely small in comparison. To be completely correct we 
would have to include these terms in the modeling problems to follow. The trouble 
is that each of the coefficients needs to be associated with a physical parameter. 
We will see that Rj 3 is physically associated with a material parameter of the porous 
medium, but it is presently unclear what the physical interpretations are for the other 
coefficients. Neglecting these terms simply leaves the door open for future modeling 

research. 
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Proceeding in a similar manner, the linearized constituent momentum transfer 



from equation (4.34) is 



J2 T a / +T j \ ={J2 ^7 + **' - £ a p aj R aj ■ v a » a . (4.38) 

\0^a J near \fi^a J eq 

The linearization constant is related to the inverse of the diffusion tensor. The lin- 



earized partial heat flux from equation 4.35 is 



^e a q a \ =K VT (4.39) 

(recalling that the partial heat flux is zero at equilibrium), and the linearization 
constant is related to the thermal conductivity. 

In each of these linearization results, the factors of volume fraction and density 
are chosen so that the linearization constants better match experimentally measured 
coefficients. The signs are chosen so that the entropy inequality is not violated. 

Several of the relationships resulting from the entropy inequality rely on proper 
definitions of the partial derivatives of the energy with respect to particular indepen- 
dent variables. The pressure is one such quantity, but there are several others that 
appear in the preceding results. For this reason, we now turn our attention to the 
exact definitions of pressure and chemical potential under our choice of independent 
variables. This will help to simplify and to attach physical meaning to the terms 
appearing in each of the linearized results. In saturated swelling porous material, 
Bennethum and Weinstein [TB] showed that there are three pressures acting on the 
system. These results are extended in the next section to media with multiple fluid 
phases. 

4.3 Pressures in Multiphase Porous Media 

We will see in this subsection that the three pressures defined in £16] can be 
extended to broader definitions in multiphase media. These definitions will help to 

simplify and attach physical meaning to the terms appearing in each of the linearized 
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results discussed in the previous subsection. We will also define several new pressures 
acting as coupling terms between the phases in multiphase media. It will be shown 
that we can return to the three pressure relationship of Bennethum and Weinstein if 
we simplify these results to a single fluid phase. 

Recall from the entropy inequality that the equilibrium pressure in multiphase 
media can be written as an accumulation of cross effects as follows: 

• -EEf^K)- (4 ' 40) 

The partial derivative is taken while holding e a , p ak , e 13 , and p^ m fixed where k = 1 : N 
and m = 1 : N, m 7^ j. Define a cross-coupling pressure as 






(4.41) 

e",p a k l£ 0, p 0n 



so that the /3-phase pressure can be simply written as the sum of these cross-coupling 
pressures 

p' s = J2p a w for /3e{l,g} and ae{l,g,s}. (4.42) 

a 

Now we derive an identity that is analogous to the three pressure relationship 
derived by Bennethum and Weinstein [TB]. To that end, consider the Helmoltz po- 
tential as a function of two sets of independent variables where there is a one-to-one 
relationship between the two sets. 

^a = ^a (^ £ a p a^ £ ^ £ /y 3 ) &nd ^a = ^a (^ p * , £ fi , p ft ) , (4.43) 

where a G {I, g} and (3 ^ a, s. The Helmholtz potential is actually a function of 
several other variables, but these are suppressed here to make the notation more 
readable. Since ip a and ip a are functions of an equivalent set of variables, the total 
differentials must be equal to each other. Setting dip a = dip a gives 
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where in each case we are taking m — 1 : N such that m ^ j, and k = 1 : N. Now 
take the partial derivative with respect to e 13 while holding e a ,p ak , and pr k fixed. In 
this case, the de a and d(e a p a: >) terms will be zero. This leaves us with: 
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Now multiply by —e a p a to get 
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Notice that the third term is p a W from equation (4.41). Define the following new 
terms: 
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(4.47) 
(4.48) 



to get the relationship 



p a w = p a w + 7r a (/3) . 



(4.49) 



Note that the new definitions only hold if e 13 ^ 0. This can be seen if one returns 
back to the Lagrange multiplier equation (at the beginning of this section) for the 
pressure. Furthermore, this relationship holds if we had taken the derivative with 
respect to e a instead of e°. 
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For completeness sake we define p 13 and n 13 so that our definitions are consistent 
with PS]: 

pP := ^p Q (« (4.50) 

a 

7r£ : =^V(«. (4.51) 

a 

With these definitions we recover the three pressure relationship derived by Ben- 
nethum and Weinstein 

pP = pP + tiP. (4.52) 

The physical meaning of p 13 is the change in energy with respect to changes in 
volume while holding mass fixed. In terms of extensive variables this is the same 
definition as pressure encountered in classical thermodynamics for a single phase. 
For this reason we call p 13 the thermodynamic pressure. The physical meaning of n 13 is 
the change in energy with respect to changes in saturation while holding the densities 
fixed. This pressure (or swelling potential as it is called in [16]) relates the deviation 
between the classical pressure, p 13 , and the thermodynamic pressure. It can be seen as 
a preferential wetting function that measures the affinity for one phase over another. 
With these physical considerations in mind we now return to the coefficient of the i? 
terms in the entropy inequality. With the present definitions, the coefficient is 



Using (4.52) this is clearly jr. Since the time rate of change of volume fraction is 
taken as a constitutive variable, the linearization result for this term can now be 
stated as 



V 



f 



n.eq. 



+ te 13 . (4.53) 



eq. 



The coefficient r arose from linearization and is formally defined as 



dp 13 



eq. 
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Equation (4.53) does little to make clear the exact meaning of this equation. The 



exact meaning will become clear in Chapter [5] under the assumption that the fluid- 
phase volume fractions are not independent. 

The definitions of the three pressures allow us to attach more physical meaning 
(and more convenient notation) to the results found when building constitutive equa- 
tions in the next sections. Before building these equations we define the upscaled 
chemical potential for a multiphase system, and after this point we will have all of 
the tools necessary to derive the new constitutive equations. 

4.4 Chemical Potential in Multiphase Porous Media 

Chemical potential is defined thermodynamically as the change in energy with 

respect to changes in the number of molecules in the system [4J [21] . This classical 
definition has the following characteristics [15] : (1) it is a scalar and measures the 
energy required to insert a particle into the system, (2) its gradient is the driving 
force for diffusive flow (Fick's law), and (3) it is constant for a single constituent in 
two phases at equilibrium. In [69], Bennethum proposed a definition for chemical 
potential in saturated porous media that satisfies all three of these criteria: 



fj?* = ^ a + p c 



dtp c 



dp c 



d {p a ip a ) 



e",p u 



dp° 



d (e a p a tfj c 



e",p° 



d {e a p a i 



(4.54) 



e",p" 



for m = 1 : N and m ^ j. In saturated media, if the changes in energy in the solid 
phase due to changes in liquid density are assumed to be zero, then the numerator of 



the right-hand side of (4.54) can be seen as the total energy in a saturated system. 



Under this assumption, the chemical potential can be rewritten as 

dijj T 



/' 



d (e a p° 



(4.55) 



e a ,p a ™ 



This indicates that in a saturated porous medium, the chemical potential of the j th 

constituent in the a— phase is the change in total energy with respect to changes 

in mass of constituent j. We now extend this definition to multiphase unsaturated 

systems. 
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Extending this idea to multiphase and multiconstituent media, we define chemical 
potential to be the change in total energy with respect to changes in mass in the 
constituent. In multiphase media we cannot make the assumption that the energy in 
one phase is not effected by changes in other phases. With this in mind, we recall that 
the total energy can be given by ipT = 'Yl a £ a p a ' l P a - Therefore, the present definition 
of chemical potential is 



fi 



Pi 



dipj 



0(eVO 



£ a ,6? ,p a k ,p/3n 



£ 



d (e a p a ip c 



d(sPph) 



(4.56) 



£ a ,eP ,p a k ,p@m 

where again, k = 1 : N and m = 1 : N, m ^ j. Notice that if dip a /dp l3j = for a ^ 



then this definition collapses to equation (4.54). Furthermore, recalling the definition 



of the Lagrange multiplier, A Qj , from (|4.15|), equation (J4.56 ) can be rewritten as 

= 1/)? + \b . 



a \ 



e a p a dip a 



e? dph 



(4.57) 



Equation (4.57) only holds for j3 = I and /3 = g. A definition of the solid phase 



chemical potential is beyond the scope of this work. 

As a result of this definition of chemical potential we observe an immediate effect 



on the rate of mass transfer terms in the entropy inequality. Using equation (4.57) 



the last three terms in the entropy inequality, (4.13), can be rewritten as 

TV 

-E4'{( A ' J +^)-( A "+^) 



i=i 



+i((^) 2 -K-) 2 ) + i((^) 2 -K^) 2 )} 
N r ii 

- E ^ ( A * + V) - ( XSJ +^ + 2 iv9,S)2 + 2 ( {V9J,9)2 ~ ( " SJ ' S)2) 

N ( 

- E* K - **> + \ M ) 2 - k ,s ) 2 ) + \ ((^) 2 - (v^f) 

- i>' |( ASj + r) - & - \ (v l >f + \ (Vo 2 - (v^) 2 ) } 

.7 = 1 ^ J 



A' 
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- r it i 

- E % ^ ~ ( A ' J + ^ + 9 (v9 ' S)2 + 9 ^'^ ~~ ^^ • (458) 

The square of the relative velocities are likely zero since these models are designed 
with creeping flow in mind. With these simplifications, the mass transfer terms from 
the entropy inequality are rewritten as 

N N N 

j=l j=l jf=l 

(4.59) 

At equilibrium we assume that the mass transfer between phases is zero. Take 
note that this is an assumption about how the constitutive variable behaves at equi- 
librium and not an assumption about the equilibrium state itself. This fine point 
is made since in several works this assumption is made as part of the definition of 
equilibrium (for example, [H]). In the author's opinion this is a subtle mistake. The 
assumption that ig = at equilibrium implies a final equilibrium relationship; the 
mass transfer between the fluid phases is proportional to the difference in chemical 
potentials 

e l i \ n , q = % U + [{P h - P 9 >) M] (ft* - //0 

= [{p h ~ P 9j ) M] (fi h - //') . (4.60) 

This helps to verify our choice of upscaled chemical potential by satisfying the third 
criteria set forth at the beginning of this subsection. Furthermore, this suggests a 
natural coupling between the liquid and gas phase mass balance equations. The mass 
transfer coefficient is chosen to have a factor of the difference in densities so as to 
better match experimental measures [7S]. Given that the units of the rate of mass 
transfer are [ML _3 t _1 ] we see that the units of the linearization constant are 
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A further verification that we have properly defined the multiphase chemical 
potential correctly can be seen through the Gibbs-Duhem relationship from thermo- 
dynamics [21]. Simply stated, the Gibbs-Duhem relationship states that the Gibbs 
potential of the a phase is the weighted sum of the chemical potentials: 

N 

r a = J2 caj v aj - ( 4 - 61 ) 

3=1 



Equation (4.61) specifies the relationship between the Gibbs potential, T a , and the 



chemical potential. Substituting (4.57) into the right-hand side of (4.61), carrying out 



the summation, and applying the definition of pressure, (4.28), gives the equation 



T a = r + - a , (4-62) 

which is the standard thermodynamic relationship between the Helmholtz potential 
and the Gibbs potential. This clearly demonstrates that the definition of multiphase 
chemical potential used here is consistent with the classical thermodynamic definition. 
At this point we turn our attention toward using the relationships derived from 
the entropy inequality to develop novel expressions for Darcy's, Fick's, and Fourier's 
laws of flow, diffusion, and heat conduction. For a concise summary of all of the 



results derived in this chapter, see Appendix D 



4.5 Derivations Constitutive Equations 

In this section we derive general forms of Darcy's, Fick's, and Fourier's laws 

based on the HMT results in the previous sections. These equations will be coupled 

with mass and energy balance equations to form a macroscale model for heat and 

moisture transport for unsaturated media. The results derived in this section extend 

the classical forms of each of these laws. These extensions suggest terms that, in the 

author's knowledge, are previously unreported. Also, we propose new forms of these 

laws in terms of the macroscale chemical potential. This suggests that the chemical 

potential is a generalized driving force for flow, diffusion, and heat transport. 

62 



4.5.1 Darcy's Law 

In 1856, Henri Darcy proposed his empirical law governing flow through saturated 
porous media |3lj . This was derived through experimentation on sand filters used to 
purify the water in the fountains of Dijon, France. In its simplest form, Darcy's law 
states that the averaged fluid flux is proportional to the gradient of hydraulic head 
(or fluid pressure) 

e l v l ' s = -kVh. (4.63) 

Under the construct of Hybrid Mixture Theory, Darcy's law is obtained by coupling 



the momentum balance equation for a fluid phase, (3.29), with the linearized constitu- 
tive equation for the momentum transfer from other phases. This has been illustrated 
by several authors (some examples include [121 HH1 IH3 EI]), and depending on the 
set of independent variables postulated for the Helmholtz potential, the momentum 
transfer term can suggest different forms of Darcy's law. 



In the present case, we recall from equation (4.36) that the linearized momentum 
transfer terms can be written as 



N 



p d P h p d P ?i ' p 



_ £ /y y ?WLvp» - etp^Vr - eV^C : V (C 8 ) , (4.64) 
F *r? dffi dJ s dC v= ; 

where we recall that R a is related to the resistivity of the medium and arose from 
the linearization process. 

Linearization of the stress-pressure relationship for the fluid phases gives an ex- 
pression for the stress near equilibrium: 

f = -p^l + ^ :g. (4.65) 
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The fourth-order tensor multiplying the rate of deformation tensor can be simplified, 
in most cases, to correspond to the viscosity of the medium (see any text on contin- 
uum mechanics). Ignoring the acceleration terms in the momentum balance equation 



(3.29), and substituting equation (4.64) for momentum transfer and (4.65) for the 



stress tensor gives the following generalization of Darcy's law: 

+ E[(^ + <)v/.]-^Egv.. 

- eVfjJvj' - £ y |^ : V (£*) + V ■ (g : g\ . (4.66) 

To arrive at this form of Darcy's law we have also assume that VC^ ~ since it 
assumed that concentration gradients in the solid phase do not affect flow. The first 
term indicates that flow is primarily due to pressure gradients, as expected. The 
eighth and ninth terms were previously reported by Weinstein in |81j . Note that the 
extra factor of e 13 on the left-hand side of the equation can be moved to the right. If 
all but the first term on the right-hand side are then ignored we arrive at the classical 
Darcy's Law 

j|/3 . ( £ V' S ) = -V/, (4.67) 

where q 13 = e^v^' s is known as the Darcy Flux. 

The linearization constant, Iff, is related to the resistivity of the porous medium, 
the inverse of which is assumed to exist, and we define K_ = (R ) ■ The tensor 
Kr is related to the hydraulic conductivity. To determine the exact meaning of the 
linearization constant we consider the units of the simplest terms: 

qP re -J^ ■ V/. 

The units of the Darcy flux are length per time [L/t] , and the units of the pressure are 

mass per length per time squared [M/(L ■ t 2 )]. This indicates that the linearization 
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constant has units [(L 3 ■ t)/M], which can be rewritten as [(L 2 )/(M/(L ■ £))]. The 
numerator of this fraction has units of permeability, k, and the denominator has units 
of dynamic viscosity, //g. This suggests that 

K P = — > ( 4 -68) 

— H 

and this relationship is confirmed in equations (11.4) and (11.5) of Pinder et al. [62J. 
The hydraulic conductivity of a porous medium is defined as 

//ore ore 
k = -^ = ^= (4.69) 

tip "P 

where v$ is the kinematic viscosity of the fluid. This indicates that Kf can also be 

defined as 

k 

K? = -it-. (4.70) 

— pPg 

It is clear from these relationships that K ^ is a function of both the type of fluid and 
the geometry of the porous medium. This coefficient "describes, in some sense, the 
ability of the porous medium to transmit fluid" [62] . In saturated porous media, the 
permeability is typically assumed only to be a function of geometry. Under Hybrid 
Mixture Theory we must note that the permeability is a function of any variable which 
is not necessary zero at equilibrium. Typically it is assumed that the permeability of 
an unsaturated medium is a function of the volume fractions [SI |62] . The tensorial 
notation may be dropped in isotropic media, but for anisotropic media it is assumed 
that the permeability may depend on the direction of flow. We will expand upon this 
idea in later chapters when building a macroscale mass balance model. 

4.5.2 Darcy's Law In Terms of Chemical Potential 



Equation (4.66) couples all of the physical processes that we wished to model at 



the outset; multiphase flow with constituents in each phase and a deformable solid. In 

order to build reasonable models based on this constitutive equation, functional forms 
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for the wetting potentials, -k^w and 7r^M , and the changes in energy with respect to 
density are needed. The solid-phase terms are likely negligible for non-deformable 
media, but dealing with the remaining terms represent a significant modeling task. 
The goal of this subsection is to greatly simplify this model while maintaining the 
physical interpretation. This is done by switching thermodynamic potentials. 

Recall from thermodynamics that the Gibbs potential, r a , can be written in 
terms of the Helmholtz potential, ip a , as 



T a = ip a + 



P 
P~ c 



(4.71) 



Taking the gradient of the Gibbs potential, expanding the resulting gradient of 
Helmholtz potential in terms of the constitutive independent variables, and multi- 
plying by —e^p 13 yields the equation: 



N „ ,a N 



E«V^. + £«V^v,-» 



3=1 



dpfa 



i=i 



dpK 



r/3^/3 



eYVT? - S -^rVp 3 - eWVT. 



P 



(4.72) 



Matching the common terms between (4.66) and (4.72), and recognizing the resulting 
chemical potential terms yields 



e /»j|.(eV'') 



pp 

7 = 1 ^ = ' 



Observe that the summation in eq. (4.73) can be simplified to 

N 



N 



J2 [a? (n* - /) Vp&] = e? (T 3 - ^) Vp? + eV £ (^ VC^) 

3=1 3=1 
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by expanding the gradient of p^ and using the Gibbs-Duhem equation (4.61). Sub 



stituting this back into (4.73) and simplifying yields the chemical potential form of 
Darcy's law: 

= _ e /y JT (C^V^) - eV^VT + e?{fg + V • (g : ^) (4.74) 

i=i ^ = ' 

Cancelling the factor of ^ from the left-hand side, rewriting the coefficient of the 
V/i^ j term, and multiplying by the inverse of Rf gives 

N 



eW 



K 



P 



J2 (p ft ! V/0 + /^VT - /<? - lv ■ (V : #) 



(4.75) 



Equation (4.75 ) states an amazing fact: the flow of phase /3 is due only to gradients 



in chemical potential, temperature, gravity, and viscous forces. The viscous forces 
are often neglected in creeping flow. This gives 



eW = -K p 



N 






(4.76) 



It should be emphasized that no additional assumptions were made to arrive at this 
equation. That is, we still assume multiphase flow with a possibly swelling solid 
phase. All of the actions, interrelations, and cross coupling effects are tied up within 
the chemical potential term. This further indicates that the chemical potential is a 
generalized force that, in effect, incorporates several driving forces. 

A final simplification is to consider a pure fluid phase where only one constituent 
is present. In this case, Darcy's law is rewritten as 



eV'' = -Kf ■ [/ Vr^ + p^rfVT - p^g] 



(4.77) 



where V 13 is the macroscale Gibbs potential. The entropy coefficient of the gradient 
of temperature poses a significant modeling issue as the entropy is not readily mea- 



surable. The fact stated by equation (4.77) is that the Darcy flux of a pure species 
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is truly controlled by gradients in temperature and Gibbs potential. This is a gener- 
alization of the classical pressure formulation that captures a wider range of physical 
effects. 

4.5.3 Fick's Law 

We now turn our attention to diffusion and Fick's law. In 1855, Adolf Fick 

published the first mathematical treatment of diffusion [37]. The empirically based 
equation simply states that the diffusive flux of a species through a mixture is propor- 
tional to the gradient in concentration of the species. This has since been generalized 
through thermodynamics and physical chemistry [211 ES] to state that the diffusive 
flux is proportional to the gradient in chemical potential of the species. In this sub- 
section we apply the Hybrid Mixture Theory construct to derive a version of Fick's 
law for multiphase porous media. It should be noted here that the classical chemical 
potential from the thermodynamic definitions of Fick's law for diffusion in a liquid 
not in a porous medium is the not the same chemical potential as that defined for the 
porous media. In mixture theory we view the porous medium as a mixture of phases 
(and species), but the classical thermodynamic definition considers one phase with a 
mixture of species. With this difference in mind, it is not immediately clear that the 
multiphase version of Fick's law will be the same. 

To derive the present version of Fick's law we first consider a linearization of the 
coefficient of the V«; aj ' a term in the entropy inequality. The gradient of the diffusive 
velocity, Vi>° J,a , is taken to be zero at equilibrium so this coefficient is zero (since 
entropy generation is minimized at equilibrium). Therefore, 

e a t a i + e a p a i (\ a J - ip a i) I + e a p a i\ aN = for all j. 



Using equation (4.32) for the definition of the Lagrange multiplier, X° N at equilibrium 



and linearizing the coefficient of Vi> aj '° about Vv aj ' a gives 

e a fi = £ a g a > : Vv a i> a + (-e Q p Qj A Qj + e a p a3 ^ aj ~ e a p a ^ a ) |, (4.78) 
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where S aj is a fourth-order tensor that arises from linearization. Now consider the 



species conservation of momentum equation (3.28). We ignore the inertial terms 
since diffusion is assumed to be slow (this is discussed in some detail in Chapter [2]). 
Now eliminate the momentum transfer terms using the linearized momentum transfer 



derived from the entropy inequality, (4.38), use (4.78) for the stress tensor, and using 



the fact that /i aj = A Qj + ip a gives a generalized form of Fick's law: 

- e a p a > ' V/x a ^ + V • ( e a S a J : Vv^A + e a p a >g. (4.79) 



The term containing the gradient of diffusive velocity is likely negligible as it is second 
order. If not, we would have to relate the fourth-order tensor, S_ a3 , with some physical 
process (similar to viscosity for fluid flow). If we neglect this term then Fick's law 
can be written as 

e«p«ilp . v a >' a = -p^V^ + p a >g. (4.80) 

Despite the novel choice of variables for this work, this form of Fick's law is identical 
to that found by Bennethum and Murad J15J and Weinstein [81| . 

The linearization coefficient in Fick's law has a similar meaning to that of the 
resistivity tensor in Darcy's law. In this case, though, we wish to associate the inverse 
of this tensor with the diffusivity tensor from classical Fick's law. Assuming that the 
inverse exists we have 

£ a p a 3v a 3 , a = _ p oyj^ . [V/i^ - g] . (4.81) 

The units of the left-hand side are [ML~H~ X \, and the left-hand side term is commonly 
known as flux. Therefore, the units of D aj is simply time [t] . Typically the diffusivity 
constant in a gas is measured as [L 2 t _1 ], so we correlate D 9j to the diffusion coefficient 
for that phase, D 9 , via the relationship 

£>* = (-=^-=) D\ (4.82) 
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where R 9j is the specific gas constant for constituent j. The units of R 9j T are [L 2 £~ 2 ] 
and the units of D 9 are [L 2 £ -1 ], hence making the units of D 9] [t\. This is consistent 
with the forms of Fick's law from thermodynamics and physical chemistry [2TJ [56] . 
Hence, the gas phase form of Fick's law is 

e?(Pv»* = - (j^f) M ■ ( V ^ J " 9] • (4-83) 

To close this subsection we finally recall from our discussion of pore-scale diffusion 
(see Chapter [2]) that the diffusive velocities are related via 

JV 

J2p aj v aj > a = 0. (4.84) 

i=i 

Multiplying by the volume fraction and recognizing the left-hand side of Fick's law 
indicates that 

£{(i£^= ' [V ^~4 = ° (485) 



near equilibrium. Equation (4.85) simply states that the gradients in chemical po- 
tential are not independent of each other. This fact will be used in future chapters 
as part of a moisture transport model. 

4.5.4 Fourier's Law 

The final result in this chapter is an extension to Fourier's Law for heat conduc- 



tion. Notice that in the chemical potential form of Darcy's Law, (4.76), there is a 
term that involves the gradient of temperature. That is, the Darcy flux is partially 
driven by a gradient in temperature. This means that Darcy flow is naturally driven 
by gradients in temperature as well as gradients in chemical potential. To properly 
handle this coupling we can either assume that the gradient of temperature is zero 
(constant temperature) or consider the energy balance equation and track tempera- 
ture as well as chemical potential. To move toward a closed system of equations, we 
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derive a version of Fourier's Law from the entropy inequality so that we have an ex- 
pression of heat flux in the energy balance equation. At the outset we first recall that 
in the entropy inequality we've assumed only one temperature for the entire porous 
medium. This implies that we've assumed that the separate phases are in thermal 
equilibrium. For this reason, we will develop an analogue to Fourier's Law that holds 
for the entire (bulk) medium. 

Following Bennethum and Cushman's work on heat transport in porous media 



[Lf] we observe that if we sum the energy equation (3.35) over a we obtain the bulk 
energy balance equation 

De 



p— t : Vv - V • q - oh = 0, 

H Dt = H H 



(4.86) 



where 



p = j>v 



ov 



J^eVV 



u a = v a - V 



t = J2 [ £ T ^ £ V« Q ® u a ] 



Er a q q i acta ol 

[s p e + e p u ■ u 

a 
a 

ph = ^e a p a h a . 



e a q a + t a -u a - p a u a [ e a + -u a ■ u c 



(4.87a) 


(4.87b) 


(4.87c) 


(4.87d) 


(4.87e) 


(4.87f) 


(4.87g) 



Given identities (a) - (g), the derivation of (4.86) follows after some significant algebra 



Define the medium velocity, v, as the weighted velocity of the medium, and the 

relative velocity, u a = v a — v, is the a— phase velocity relative to the medium. Note 

that v a,s = v a — v s — v+v = u a — u s = u a,s . In the case where the velocity of the solid 

phase relative to the medium is zero (u s = 0) we immediately see that u a = v a,s . This 
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assumption along with equation (4.87f) indicates that there is naturally a coupling 



between the relative velocities, v a,s , and the total heat flux, q. 



Using the near equilibrium result, ^2 a £ a q a = K_ • VT (equation (4.39)), we can 
write the total heat flux as 

1 



q = K-VT + J2 



/-'» ,„'!•■ _ ft 01 ** 01 ' 3 I ■ " 



t'v' 



e + r 



a,s „,a,s 



(U 



If we were to (wrongly) neglect all of the terms in the summation we would arrive 
at Fourier's Law for heat conduction. The trouble here is that the terms in the 
summation are not negligible, and therefore the total heat flux in a porous medium 
must be a function of the gradient of temperature, the relative velocities, the stress 
in the fluid phases, and the internal energy. 



Since v s ' s = v s — v s = the right-hand side of equation (4.88) is only a function 
of the fluid velocities relative to the solid phase. Neglecting viscous terms we recall 
that the fluid-phase stress tensors can be rewritten as t a = — p a J. Neglecting the 
second-order term, v a,s ■ v a,s , the total heat flux is now written as 

q = K ■ VT - J2 l(P a + P° ea ) ^"'1 • ( 4 - 89 ) 

a=l,g 

At this point we replace the internal energy term with Gibbs energy in hopes of 
deriving an extended Fourier's Law in terms of the chemical potential. Recall from 
thermodynamics that the Gibbs potential and internal energy are related through 



V 
T a + Tri a -^ 

P c 



(4.90) 



Therefore, the total heat flux can be written in terms of the Gibbs potential as 



q = K ■ VT - Y^ [P a ( F " + T V a ) v c 

a=l,g 



(4.9i; 



Using the Gibbs-Duhem relationship, (4.61), this can be rewritten in terms of the 
chemical potential as 



K ■ VT - Y 

a=l,g 



N 






(4.92) 
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The trouble with both (4.91) and (4.92) is that they both rely on measurements of 



entropy. One way to work around this issue is to assume that the entropy is only a 
function of temperature, and then to recall that the specific heat is defined as 



T 



drj a drj 



T- 



dT 



Solving this separable ordinary differential equation (under the assumption that the 
variation of specific heat with temperature negligible) gives 



V a (T) = c p ln(^)+r ] % 



(4.93) 



where To is a reference temperature, and r/g is a reference entropy. While this is only 
an approximation it does allow us to move forward without direct measurements of 
entropy. 



The extended Fourier's Law (4.91) presented here frames the equations presented 



in [14] in terms of the Gibbs potential. This will allow for easier coupling with 
the chemical potential forms of Fick's and Darcy's Laws presented in the previous 



subsections. The caveat is that the equation for total energy balance, (4.86), is not 



particularly useful since we do not have constitutive relations for the total stress and 



total energy. For that reason, we will not use equation (4.91) or (4.92) for Fourier's 



law in the energy equation. Instead we will use the linearized partial heat flux and 
the constitutive relations for the phase stresses and relative velocities to derive a 
generalized heat equation. 

4.6 Conclusion 

In this chapter we have shown that a novel and judicious choice of independent 

variables for the Helmholtz Free Energy can be used to derive forms of Darcy's, 

Fick's, and Fourier's Laws for multiphase porous media. These equations are similar 

to those found in [TU HH [151 El]- Each equation can be written with an eye toward 

the macroscale chemical potential, and in each case the chemical potential form is 
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more mathematically appealing in the sense that there are fewer terms and many 
of the physical processes are manifested in the chemical potentials. This illustrates 
the usefulness of the chemical potential as a modeling tool. Furthermore, since the 
chemical potential appears naturally in each of these equations we have set the stage 
for a more natural method of coupling the fluid flow, diffusion, and heat transport. In 
Chapters [5] and [7] we will couple these equations with the upscaled mass, momentum, 
and energy balance equations to yield a system of equations that will govern total 
moisture transport and heat flux in unsaturated porous media. 
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5. Coupled Heat and Moisture Transport Model 

To form governing equations for heat and moisture transport in porous media we 
pair the constitutive equations derived in Chapter [4] with upscaled mass, momentum, 
and energy balance equations derived in Chapter [3] There are several existing mod- 
els for each physical process of interest (fluid flow, diffusion, and heat transport) and 
recent research indicates a need to understand the fully coupled system of equations 
as it relates to moisture transport, evaporation, heat transport, and other physical 
phenomena. In this chapter we derive a model for coupled heat and moisture trans- 
port using Hybrid Mixture Theory and knowledge of pore-scale effects. To begin this 
modeling task we first investigate the classical models used within the past century in 
Section 5A_ In Sections 5J2 and 5J3 we pair our constitutive equations from Chapter [4] 



with upscaled balance laws from Chapter |3j perform a dimensional analysis, and dis- 
cuss forms of the linearization coefficients arising from HMT. This is done in an effort 
to generate a closed system of governing equations. Several simplifying assumptions 



are made to close the system in Section 5.4 The solution(s) to the closed system will 
be discussed in Chapter [7j 

5.1 Introduction and Historical Work 

To give the reader a better understanding of the work from the past century, we 
present three classical models here with some discussion on their advantages and dis- 
advantages. First we discuss Richards' equation for unsaturated fluid flow in Section 



5.1.1, second we discuss Phillip and De Vries enhanced diffusion model in Section 



5.1.2, and lastly we discuss De Vries' heat transport model in Section 5.1.3 



5.1.1 Richards' Equation for Fluid Flow 

The classical equation for fluid flow in unsaturated media is known as Richards' 
equation (also called the saturation equation). This equation was first derived in 1931 
by L.A. Richards at Cornell University [65]. It takes a postulated form of the mass 
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balance equation (similar to equation (3.25)) and replaces the flux term with Darcy's 



law. The gradient of pressure is rewritten in terms of pressure head (h = p/(pg)), 
and then a constitutive relation is assumed for the pressure head as a function of 
saturation (or volume fraction). Another constitutive relation relating the relative 
permeability of the medium to saturation is assumed. There are several versions of 
the constitutive relations, but one of the more popular in recent research are those 
of van Genuchten [791 162] . Another more recently investigated relationship is the 
Fayer-Simmons model J36J EHJ [78], which is an extension of the van Genuchten model 
to cover the case of very low saturations. 

The result of the assumption and substitutions in the mass balance equation is a 
nonlinear diffusion equation where the primary unknown is the percent saturation of 
the medium 

— = V-[D(S)VS-K(S)z], (5.1) 

where K(S) is the hydraulic conductivity function and D(S) is the product of K(S) 
and the derivative of capillary pressure with saturation. Recall that saturation is 
defined as 



e l e l e l volume of liquid 

1 — e s e 9 + e l e volume of pore space 



S = -?— = -^—f = -= VU1UlliC Ui 114U1U (5.2) 



and is understood as the volume of liquid per volume of pore space. 

This model has been effectively used for several decades, but there are a few dis- 
advantages of note. First of all, this equation does not allow for phase change between 
the liquid and gas. The original model was proposed for systems with immiscible flu- 
ids, where phase changes likely don't occur, but it is also used for unsaturated soils 
where phase change is possible and air is always availabe.. A second disadvantage is 
that humidity and temperature gradients are not considered. A third disadvantage 
is that the pressure head - saturation curve is hysteretic (depends on the history of 
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flow). The constitutive laws for pressure head don't account for this hysteretic behav- 
ior directly. Instead, it is often assumed that fitting parameters change with changing 
direction of flow. This leads to the final disadvantage: the use of the van Genuchten 
capillary pressure - saturation relation. This is a widely used relationship, but re- 
lies heavily on two fitting parameters. The measurement of these fitting parameters 
is difficult, and they are typically found by fitting numerical solutions of Richards' 
equation to experimental data. 

Several extensions and modifications to Richards' equation have been made re- 
cently, the most notable of which is that of Hassanizadeh et al. [191 Wf\ . In these 
papers, they propose a dynamic relationship between capillary pressure and satura- 
tion based on Hybrid Mixture Theory with interfaces. They also propose that the 
hysteretic effect observed in the capillary pressure - saturation curves is due to the 
(postulated) fact that the capillary pressure, saturation, and interfacial area density, 
e lg , form a unique surface. This partially explains hysteretic effects by seeing them 
as a projection of this surface onto the capillary pressures - saturation plane in the 
p c — S — e l9 space. This model is gaining in popularity, but is far from widespread 
acceptance. Some of the relevant publications are [Ml HH1 [501 SHI H?l [58] . 

In the present chapter we present a modification to the Richards' equation that 
incorporates the dynamic capillary pressure relationship of Hassanizadeh et al. The 
major differences between the present derivations and their work are: (1) modeling 
in terms of chemical potential, (2) allowing for phase transition, and (3) allowing for 
humidity and temperature gradients. Our present modeling effort will account for all 
of these effects, and hence, constitutes a generalization of the existing model. 

5.1.2 Phillip and De Vries' Diffusion Model 

In 1957, Phillip and de Vries published their comprehensive work on diffusion of 
water vapor in porous media |611 . In their model they postulate an enhanced Fick's 
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law, 

q*> = - p 9Dt]VC 9 \ (5.3) 

where q 9v is the water vapor flux and 77 is an enhancement factor that is a function of 
the toruosity, volume fraction of air, and a "mass-flow factor" . The mass-flow factor 
is then postulated as a function of pore-scale gradients in saturation and temperature. 
This model has successfully been applied to several diffusion and evaporation prob- 
lems (e.g. [IB])) but the trouble is that the exact form of the enhancement is based on 
empirical evidence. Furthermore, this model has come under recent scrutiny due to 
the fact that the proposed factors affecting 77 are pore-scale effects and are therefore 
difficult to accurately measure [251IZI1E21EQ1II31EI1E81EQ]. Many of these works 
use x-ray tomography to attempt to measure these pore-scale effects directly. 

In the work by Cass et al. [21], an empirical form of the enhancement factor was 
proposed. In this work, a fitting parameter is used in the enhancement factor to 
arrive at good agreement with experimental data. This model has been used in more 
recent works (e.g. [68| [78]) in conjunction with a mass balance equation for the water 
vapor in the gas phase. The resulting model is a nonlinear diffusion equation for 
concentration of water vapor that deviates from the more classical de Vries model. 
Aside from the empirical fitting parameter, the mass transfer between phases also 
relies on a fitting parameter and an empirically-derived functional form. 

In the present chapter we build a model for diffusion based on using the chemi- 
cal potential as a primary unknown and the Hybrid Mixture Theory construct. The 
enhanced diffusion is not incorporated into these models, and the mass transfer is 
modeled by the difference in chemical potentials; a more physically natural formula- 
tion. A comparison will be made to the model of Cass et al. 

5.1.3 De Vries' Heat Transport Model 

In 1958, de Vries published a second paper coupling heat and moisture transport 

in porous media [32]. In this research, he proposed an extended heat transport 
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model for porous media that is still used today. Neither his diffusion nor his heat 
transport model were thermodynamically derived. Instead, he began each derivation 
with a postulation of the forms of diffusive and heat flux. For the heat transport 
equation he included terms similar to the classical Fourier's law, but also proposed 
that heat transport was due to advective transport in the fluid phases. This model 
is still popularly used today to couple heat and mass transport in unsaturated media 
[5j [77J [78]. That being said, the effects included in this equations are based solely on 
de Vries' supposition of the factors affecting heat flow. 

In 1999 Bennethum and Cushman published (to the author's knowledge) the first 
work using Hybrid Mixture Theory to derive an extended de Vries model for heat 
transport in swelling saturated porous media [14J. In the present chapter we take a 
similar approach using HMT to derive a thermodynamically consistent model for heat 
transport in non-swelling unsaturated media. This is done with an eye toward using 
gradients in temperature as the thermal diffusion process and the chemical potential 
to describe the secondary processes such as advection. 

5.2 Assumptions 

In this section we state the baseline assumptions that will be used throughout the 
remainder of this work. These assumptions are meant to make minimal limitations 
on the applicability of the resulting models, but at the same time they are meant to 
keep the mathematics tractable. Possible relaxations to these assumptions (and the 
source of possible avenues of future research) will be stated as they are encountered. 

The simple set of baseline assumptions are as follows: 

Assumption #1: The solid phase is rigid, incompressible, and inert. 
Assumption 7^2: The liquid and gas phases are each made up of N constituents. 
Assumption 77=3: No chemical reactions take place in any of the phases. 
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The first assumption is the most restrictive. Mathematically it corresponds to 
setting the Lagrangian derivatives of both density and volume fraction for the solid 
phase to zero. Assuming that the solid is inert simply means that no mass will 
precipitate onto, or dissolve away from, the solid phase. With these assumptions, the 



solid phase mass balance equation (from equation (3.25)) becomes 



V-v s = 0. (5.4) 

If a deformable solid is considered where the solid-phase volume fraction can change, 
then this assumption would need to be relaxed. One particular relaxation of this 
assumption is to allow for incompressibility and inertness of the solid phase but relax 
the rigidity assumption. Under this relaxation, the solid phase mass balance equation 
becomes 

D s e s 

' - - £ s V--u s = 0. (5.5) 



Dt 

A consequence of fixing the solid phase volume is that e l + e 9 = 1 — e s := e, 
where e is known as the porosity of the porous medium. A further consequence is 
that the liquid and gas phase volume fractions are no longer independent of each 
other. Note that we could have made this assumption up front and exploited the 
entropy inequality with this assumption (this is done in jHJ H5] for a different set 
of independent variables), but proceeding in this order allows us to return to the 
present entropy inequality results and consider a deformable solid in the future. Since 
the fluid-phase volume fractions are no longer independent we can replace them by 
saturation as defined by 

e l e l 

S = -^— = -. (5.6) 

e l + e9 e y ' 

This implies that the volume fractions are related via e l = eS and e 9 = e(l — S). 

Assumption #2 is a byproduct of the principle of equipresence and will be relaxed 

later for simplicity. In the most general sense, this assumption states that every 
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species that exists in one fluid phase also exists in the other. In reality this is likely 
not true. For example, if a constituent is present in the liquid phase it is possible that 
evaporated particles of the constituent are not be present in the gas phase. Another 
example would be if we were to extend this model to an oil-water system. The two 
fluids in this case are immiscible and it is unlikely that every species in the water 
phase is present in the oil phase (and visa versa). We take this into account by 
setting the appropriate concentrations to zero after the constitutive equations have 
been derived. 

Assumption ^3 indicates that the rate of mass exchange due to chemical reac- 
tions, f Qj , is zero for all phases. The consequence of this is that the rate of mass 
generation of a constituent in a phase only occurs between two phases. This is true 
for some porous media, but chemical reactions can occur in some specific cases such as 
remediation problems. Under this assumption these cases are henceforth eliminated 
from the discussion. 

Other simplifying assumptions exist for many media, but the three presented 
herein constitute a set that leads to several mathematical simplifications with as few 
physical restrictions as possible. 

5.3 Derivation of Heat and Moisture Transport Model 

In the remainder of this chapter we focus on using the results from Chapters [3] 
and |4j along with the assumptions from Section 5.2 , to derive a closed system of 



equations for heat and mass transport in unsaturated porous media. This will be 
done with an eye toward using the chemical potential as the driving force for these 
processes. We will show that under certain additional simplifying assumptions that 
a closed system can be derived. 

5.3.1 Mass Balance Equations 

We first build generalized mass balance equations in terms of the chemical po- 
tential under assumptions #1 - #3. Recall from Chapter [3] that the mass balance 
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equation for the j th constituent in the a— phase is (from equation (3.21)) 

{ ^ t P ' + e a p a *V ■ v a > = J2 e? + r a '- (5.7) 



The last term can be dropped under assumption #3 in Section |5.2[ Because of the 

form of the constitutive equation, and to adhere to the principle of frame invariance, 

it is convenient to rewrite this equation relative to the solid phase. To do so we recall 

the identities 

D a i(-) D a (-) 

_i) = _A) + ^. V (.) (5.8a) 

D a (-) D s (-) 

_M = _M + ^. V{ . ) ( 5. 8b ) 

and expand the Lagrangian time derivatives accordingly to get 

DS ( £ *P a ^ + v <*i,* . V ( £ a p a i) + v a ' s ■ V (e a p a J) + s a p a ^V ■ v a > = ^ eJ J . (5.9) 

Taking the definition of the Lagrangian time derivative, 

Dt dt + U ' 

adding and subtracting e a p aj 'V ■ v a , and subtracting e a p°' 3 'V ■ v s = gives 

r\ / Qi fy ■ \ 

[£ P J) + V ■ (eV'« ai,a ) + V ■ (e a p a iv a > s ) = J2 *p- ( 5 - 10 ) 

dt fra 

Notice the use of Assumption #1 in the last step, and observe that if Assumption 
#1 is relaxed then the mass balance equation would involve a time derivative of the 
solid-phase volume fraction (at least). 



Equation (5.10) is the general mass balance equation for both of the fluid phases. 
Notice that we are not replacing the volume fractions with saturation here since we 
don't know if a is the liquid or gas phase. Substituting Fick's law for the diffusive 
flux and Darcy's law for the Darcy flux gives the chemical potential form of the full 
mass balance equation for species j in phase a: 

{ / t } - V • {fg» ■ [Vp a > - g}} 
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V • < p a iK c 



N 



J> a *v 



/' 



"fc x 



p a rj a VT - p a g 



fc=i 






(5-11) 



N 



2(p afc V^) + pYVT-A 



-/3- 



(5.12) 



It should be noted here that the Eulerian and Lagrangian time derivatives are equal 
under the assumption that the solid-phase velocity is zero (Assumption #1). Also 
note that if we sum over all constituents then we arrive at the mass balance equation 
for the phase (where we have used X) 7 -=i p a Jv a J ,a = 0) 

dt I = 

I U=i 

The chemical potential form of the mass balance equation is only one form. We 
could have used the pressure formulation for Darcy's law and arrived at a pressure - 
chemical potential form of the mass balance equation. 

The rate of mass transfer term on the right-hand side of the mass balance equation 
can be rewritten in terms of a linearized result from the entropy inequality. Recall 



E' 



from equation (4.60) that the mass transfer term can be written as 



[(p a > -ffi)M] (ffi ' -//■>'), 



(5.13) 



where the coefficient (p aj — p^ J ) is chosen to be consistent with equation (9) of [78] . 
Also recall that since the interface is assumed to contain no mass we must have that 
the rate of mass gained from the j3 phase to the j th species in the a phase must be 
equal to the rate of mass lost from the a phase to the j th species of the j3 phase: 
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If the chemical potential of the liquid phase is larger than the chemical potential of 
the water vapor then mass will transfer from liquid to gas and e l < 0. Similarly, if 
the chemical potential of the liquid phase is smaller than that of the water vapor then 
mass will transfer from gas to liquid and e l > 0. Recall from the discussion adjacent 



to equation (4.60) that the units of M are the reciprocal of flux. 



S3 



There are clearly more unknowns than equations in the 2N fluid equations since 
we must account for the densities, temperature, volume fractions, and entropies as 
well as the chemical potentials. Certain sets of simplifying assumptions can be used 
to reduce the number of unknowns (e.g. incompressibility of a fluid phase). These 



will be discussed in Section [5\4| Instead of making these assumptions up front we now 
turn our attention to deriving a generalized energy balance equation to account for 
the temperature. This will give one more equation but will add no more unknowns 
to the system of equations. 

5.3.2 Energy Balance Equation 

As another step toward developing a closed system of equation equations for 

heat and moisture transport we next examine the energy balance equation. This will 
give an equation in terms of temperature, chemical potentials, saturation (volume 
fractions), entropy, and densities; increasing the equation count but not increasing 
the variable count. Since we assumed at the outset that all of the phases are in thermal 
equilibrium we will only have one equation for energy balance. This will be derived 
by considering the sum of each of the phase energy balance equations. Counter- 



intuitively, we will not use the form of Fourier's Law (equation (4.87f) or (4.92)) 
derived for the total heat flux since the energy equation derived in that section is 
more cumbersome to work with than the individual phase energy equations. Instead 
we will use the partial heat flux for each phase as derived from linearization about 



equilibrium (4.39). 



From equation (3.35), the volume averaged energy balance equation is 



e a p a ^- - e a f : # - V • (e a q a ) + e a p a h a = ^ Q%. (5.14) 

Using the identity Dt = —J£- + v a,s • V(-) and using dot notation for material time 
derivatives allows us to rewrite the energy equation as 

e a p a e a + e a p a v a ' s ■ Ve° - e a f : £ - V ■ (e a q a ) + e a p a h a = ^ Q% (5.15) 

feat 
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The trouble with (5.15) is that the first and second terms contain the interal energy 



density, e a . To tie this equation back to the HMT framework we've used throughout 
(and to give the equation a more natural set of dependent variables) we perform a 
Legendre transformation to change the energy term into the Helmholtz potential via 
the thermodynamic identity e a = ip a + Ti] a . The energy equation is now written as 



Dt 



+ e a p a v a ' s ■ Vip a + e a p a Tr] a + e a p a Tv a > s ■ Vr] 



+ e a p a r] a f + e a p a r] a v a ' s ■ VT - e a f : £ - V ■ (e a q a ) + e a p a h a . (5.16) 

Next we seek to remove the Helmholtz potential and entropy terms from the 
energy equation. To do this we recall that the Helmholtz potential is a function of 
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all of the variables listed in (4.6). Under the assumptions listed in Section 5.2 



drop the solid phase terms from this list. Furthermore, we know that under these 
conditions the volume fractions are not independent so we could replace both e l and 
e 9 by saturation, S. This is not done (yet) as the entropy inequality was exploited 
while assuming that they are independent. The switch can be made at any point 
later. Therefore, under the present assumptions, 

^ = ^ { £ \ e \ p h 5 p9i ? T ) forj = i : N. 



Entropy, r) a , is assumed to be a function of the same set of variables (since rf 1 = 
—dijj a /dT). Using the chain rule to expand all of the derivatives of tp a and rj a in 



equation (5.16) we arrive at an expanded form of the energy equation: 
'dip c 
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(5.17) 



From the entropy inequality we know that the temperature and entropy are con- 
jugate variables. For this reason we can cancel these terms from the T and VT 
coefficients. 



Equation (5.17) is an expression of energy balance for phase a, but since we are 



working under the assumption that the phases are in thermal equilibrium we now 
sum over all of the phases to form one energy balance equation for the entire porous 
medium. The sum is: 
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(5.19) 
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The T term can be rewritten as Y^ Q {e a p a T -^r} T = pc p T, and in doing so we 
implicitly define the volumetric heat capacity of the entire medium: 






Next we recall from equation (4.39) that the partial heat flux can be written as 



Y2a { ea Q a } — K. " ^T (more will be said about the functional form of K^ in future 
sections). The heat source term can be rewritten as Y^ Q {s a p a h a } = ph, where h is 
any internal source or sink of heat on the entire medium (i.e. heat sources that are 
not boundary conditions). 

Notice that several of the gradient terms are the same as those in the linearized 



constitutive equation for the momentum transfer, (4.31) and (4.36). Replacing these 



terms with the remainder of the momentum balance terms and simplifying gives 

E | E Qp } = p°p f - v • (£ • VT ) + p h - E i £ T ■ i a } 
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Equation (5.20 ) expresses the energy balance for the bulk porous medium. Several 



of the terms can be simplified at this point. Toward this goal, we will 



1. derive a relation for the energy transfer terms: Y2 a ^2s^a Q 



2. rewrite the stress term, E^ Q £ Q £ a '■ $*i using constitutive relationships for t a 

3. rewrite the momentum transfer terms, T , using the linearized momentum 



transfer from the entropy inequality, (4.36) 



4. rewrite the advective terms, v^ ,s , using Darcy's law, and 

5. relate the changes in entropy, #y, to material coefficients. 

The first two of these are discussed in the following two subsections. The third and 
fourth come as a consequence of the first two, and the fifth will be discussed under 
proper simplifications in future sections. 

5.3.2.1 Energy Transfer in the Total Energy Equation 

Consider the energy transfer and stress terms: Qp , Q a: > ,&n.d t a . From equations 



(|3.36a|) and (|3.36b|) we recall that the restrictions on the interface are 

= Va, 



N 

E 



Q a > + l a ° ■ v a *< a + r a > ( e aj + -v a *' a ■ v a >' c 



EE 

a Pj^a L 



Q" J + 2V • v a * + e" J e a > + -v a J ■ v 



P 



J n> a i I Z a J I n a i _L n, a i », a i 



(5.21a) 
j = 1: N. (5.21b) 



We also note the identity 

N r 



3=1 L ^ ' 



(5.22) 



(see Appendix A. 2 of [HI])- With these three identities, the sum of the energy transfer 
terms can be written as 



E E Qfi = E E E fe + f 7 ■ ^ + % U^ + \ v ^ ■ V *A 
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(5.23) 



Next we examine the momentum transfer term appearing in equation (5.23). 



Recall from equation (3.33) that 



N 
3=1 



' + ^v a ^ a 



(5.24) 



Rearranging this identity and multiplying by the a— phase velocity we see that 



N N 



3=1 3=1 



(5.25) 



Substituting (5.25) into (5.23), simplifying, and neglecting the second-order terms in 



velocity we see that 

EE{«s} = -E{^- M }-E{^-""}-4(»'-^)- (5-26) 

a p^a fcl (3^g 

Notice from this simplified version that we have eliminated the energy transfer in 
favor of the mass and momentum transfer terms after summing over a (and neglecting 
second-order effects). 



5.3.2.2 Stress in the Total Energy Equation 



We next derive the proper form of the stress term in equation (5.20). The 
a— phase stress near equilibrium is given by t a = —p a I_ + v a : d a from the lin- 
earization of the fluid phase stress tensors about equilibrium. For the solid phase 
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stress tensor, on the other hand, we will not use constitutive relations for t s but keep 



in mind that it is the sum of effective and hydrating stresses (see equation (4.23)). 
Therefore, 

J2 Kf : €) = J2 y \-p a l + ^ ■€) '■£(+ £ % s : i s 

= - J2 { £a p a l ■ €) + Yl IV '■ € : i a \ + £S i s : i s ( 5 - 27 ) 

a=l,g a=l,g ^~ 

The second term is likely negligible as the viscous terms typically play little role in 
creeping flow. This means that Yl a {^ a t a '■ d a } can be approximated by 

J2 {e°t a '■£} = ~Y1 { £a P a l '■£}+ £ % S : i S ( 5 - 28 ) 

a ot=l,g 

Using indicial notation we note that for the fluid phases, J : d a — I_ : (Vv a ) sym = 
SijVfi = vf t = V ■ v a , and therefore the stress tensor terms can be simplified to 

J2 {>t '■£} = -J2 { £a p° v ■ v °} + £S t '■ i s ( 5 - 29 ) 

a a=l,g 

The solid phase rate-of-deformation tensor is related to the strain rate of the solid 
phase. Assuming that the strain rate is zero (for a rigid and incompressible solid), 



we can neglect this term. This implies that the stress tensor term in (5.20) can be 
approximated by 

J2 Kf : £} = ~ J2 ( £ VV ■ v a } . 



Using Assumption # 1 from Section 5.2 for the divergence of the solid-phase velocity 



(V ■ v s = 0), we finally conclude that the stress term in (5.20) can be simplified to 
J2 {^T : £} = -eV V • v l ' s - eyV • v 9 > s . (5.30) 

a 

Not surprisingly, this states that the stress is related to the fluid pressures. 
5.3.2.3 Total Energy Balance Equation 
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In this subsection we use equations (5.26) and (5.30) to simplify the energy bal 



ance equation, (5.20). Substituting these into (5.20) and canceling the momentum 



transfer terms gives 

= pCp f -V-(K- VT) + ph + eVV • v l ' s + eVV • v 9 ' s + e l g (e l - e 9 ) 

- + E l" a p° 



+EWE 



cfe* de l 



de9 de9 



i=i 



dip a dr\ a 
dp l i dp 1 * 

dip a „ drf 



dp 9 * dp 9 i 



V'j 



E 

P=l,9 



N 



3=1 

drf 



E «v 



H-eV 



T 



<9T 

AT 



+E 



VT + 



T 



T 



dpv _ 



dip a 
dp 13 -* 

dr] 13 
W 

Vp^ 



+ e^p^Tp^- I Vpft 



W + 



T 



0p£ 



Ve fi 



V 



A* 



(5.31) 



Next we discuss the e a p a 'V ■ v a ' s and p a V (e a ) ■ v a,s terms. Using the product 
rule it is clear that the sum of these two terms gives p a V ■ (e a v a,s ). A choice is made 
here to remove these terms in lieu of mass transfer terms. To do so, we recall from 
the mass balance equation that 

D s (e a p° 



Dt 



+ V ■ (e a p a v a ' s ) = J2 e 



/3' 



P*a 



and solve for V • (e a v a,s ): 



p a V ■ (e a v c 



-p a e a - e a p a - e a v a > s ■ Vp Q + e%. 



We have dropped the summation on the mass transfer term since we are assuming 
that the solid phase is inert and that there are only two fluid phases. Multiplying by 
(p a /p a ) gives an expression for p a V • (e a v a,s ): 



p a V ■ (e a v c 



-p a e a 



e a p a 



a a* 



P" 



Vp° + ( j£) e% (5.32) 
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Substituting this into the energy equation gives 



=pc p f - V • (K ■ VT) + ph+ (( - x + e l 



- + e 9 



y + E i £apa 
+ \ -p g + E \ £a p a 

N 

+E 



de l 



de l 



dip a dr] 



3=1 

N 

+ E 

3=1 

+ E 

0=1,9 



£ lpl 



de9 de9 

dtp 



e g 



e 9p9 



' N 

E 

.3=1 



E < £a p a 



E i £a p a 



Q a 

\-T— — 

dp l i dp l i 

dip a i drf 



T 



eV 



dr] 13 

N 



dp 9 i dp g j 

dip a 



+E 

3=1 



VT + 
drj f 



E^v 

di] f 



+ T 



p 9j 



drf 



T 



T 



dpv 



de l 
Vp^ 



dpPj dpPj 
dr] 13 



Vff* 



Ve' + 



v 



T 



P,s 



de9 



Ve £ 



(5.33) 



There are several more simplifications that can be made. To help with these 
simplifications recall the following definitions for enthalpy, pressure, wetting potential, 
chemical potential, and entropy respectively: 



H c 



jt 



7T 



p u 



+ e c 



N 



EE 

a j=l 



e a p a p^ dip c 



E ( £a p a 



£ dpfr 

dip 



£ a p a Q^a ' 



^+E(^) 



dip a 
~dT 



With these identities in mind we make the following four simplifications: 



(5.34a) 
(5.34b) 
(5.34c) 
(5.34d) 
(5.34e) 
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1. coefficient of the mass transfer term: 



P , J 



V 



s , + e ' ) -\j 9 +e9 ))^= ( Rl ~ H °) 4 : = H 



Recalling that e l is the rate of mass transfer between the fluid phases, L is 
understood as the latent heat of evaporation since this represents the heat lost 
or gained due to phase exchanged between the fluids. This is consistent with 
the chemist's definition of latent heat as the change in enthalpy. 



2. coefficient of the time rates of change of volume fractions: 



— jf 



£ 



e a p a 



dip a ^drf 



fi + ^ + T 



-p 



-f - T 



~dT 



W 



where we recall that p 1 is thermodynamic pressure as defined in Chapter kl 



p P =p -P + 7r /?_ 



At this point we can exchange the time rates of change of volume fractions for 
time rates of change of saturation. That is, recall e l = eS and e 9 = —eS. The 
sum of the two associated terms is 



(p 9 -p l )+T 



dn 9 dn l 



dT dT 



eS. 



From the near equilibrium results from the entropy inequality we now recall 



(from equation (4.53)) that 



V 



p. 



T£ h 



n.eq. 



eq. 



(5.35) 
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Therefore, the S term becomes 



]f 



P 



cq. 



cq. 







+ 2reS + T—(^-7i l ) 



dT 



eS. 



(5.36) 



The first set of parenthesis in (5.36) (approximately) represents the capillary 
pressure as measured at equilibrium, 



P l 



<-'<!■ 



cq. 



Pc 



This will be discussed in more detail in Section 15.4.1.11 The middle term in 



(5.36) is an effect of the dynamic pressure-saturation relationship (equation 



(4.53)). The temperature derivative can be interpreted as the effect of temper- 



ature on the relative wetting potential. That is, how much does temperature 
affect the relative affinity for one phase over the other. It is likely that a con- 
stitutive equation is needed for this relationship. 

3. coefficient of time rates of change of densities: 

We wish to rewrite these coefficients in terms of enthalpy and chemical potential 
since it provides a mathematically simpler expression. 



N 

E 



£ /y 



£ < £a p a 



dip a drf 



P" 



N 

E 



: '" ; V^(/i_/)_ £ ^^_/) 



P 



P 



pP 



N 






N r 



ft _. T d J^ \ a -', 



dT 



fi 



dT 



P 



dT 



P 



where we recall that H@ is the enthalpy of phase f3. 



4. coefficient of relative velocity: 

For this coefficient we again use the definitions of pressure, chemical potential, 
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and entropy. We also rely on the Gibbs-Duhem relationship (4.61). 



N 

E 

3=1 



P 



,P 



E < £a p a 



+ eV 



T 



N 
3=1 



e p p^ 



VT + 
' T drf_ 

N 



T 



dip a drf 

drj 



VpP- 



de l 
Vp^ 



Ve' + 



P' 



T 



drj 



Ve £ 



Vp? + J] [e (^ ~ 1> p ) Vp ft ] + eV< VT 



3=1 

. d \d^„ , dip? 



dT I de l 

N 



de9 






J 
N 



-e^VpP + E [eV'V/'] + eV^VT 



3=1 



+T±L(?.+ir)+w + Y;£ 



a j=l 






JV 



3=1 

iV 



3=1 







N 



+ T " J (e^) 2 ^ ■ ^> s + p^Ve' 3 - e^^Vp' 3 + ^ [eV V/'] 



i=i 



A? 



-e^Vp' 3 + Y, [eV ft Vp^] + eV^VT 



3=1 



+ T^l (e^g ■ v^ - s^Vp? + (fy V (,V) + £ [^ ft V/'] 



Since this coefficient is contracted with the relative velocity, ir ,s , we can likely 

neglect the relative velocity term in the temperature derivative as it will result 
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in second-order effects. This simplifies the coefficient of the relative velocity to 



N 



+ T 



9 ' e^Vp^ + ^[eV J Vp ft ] + (^)v(eV) 



dT 



(5.37) 



After these four simplifications and rearrangements, equation (5.33) is now rewrit- 
ten as 



= pc p f - V • (K ■ VT) + ph + Le\ 



+ 



(V -pi )+ 2r ^+ T ^( 7r9 - 7r 

V eg. leg./ CJ-* 



£5 



JV 



- eSff'p' + e<S ^ 



i=i 



^ - T 9 * * ' 



dT 



N 



e(l-S)H 9 p 9 + e{l-S)J2 



i=i 



P 



tfi _ T ^l l ,■//, 



dT 



P* 






JV 



.^r^V/ + ]T [eV' V/'] + eV^VT 



i=i 



+T^ |_ e ^V/ + £ [eV^Vpft] + (^) V (eV) 



v 



P,s 



(5.38) 



Equation (5.38) depends on temperature, wetting potentials, enthalpies, chemical 



potentials, Gibbs potentials, saturation, densities, pressures, and relative velocities. 
Since the Gibbs potentials are functions of densities and chemical potentials this 
does not add more unknowns to the system of equations. The pressures and relative 
velocities can be paired with forms of Darcy's law, and constitutive equations are 
needed for the enthalpies and wetting potentials. We now turn our attention to the 
coupling of the fluid-phase mass balance equations and the present energy equation. 

5.4 Simplifying Assumptions — A Closed System 
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A host of simplifying assumptions can be made on the system consisting of equa- 



tions (5.11) (for a = l,g) and (5.38). These are made to reduce the number of 
unknowns and equations to a count that is more easily handled by numerical solvers. 
This is also done to avoid having to model any secondary (possibly second-order) 
physical processes (examples of which include very slow processes such as those on 
the order of (v l,s ) 2 or (v 01 ^ ) 2 ). These assumptions are in addition to Assumptions 



#1 - #3 made in Section 5.2 



Assumption tt=4: Assume that the liquid phase is composed of a pure fluid with 
no additional species. Strictly speaking this is not realistic since the water in 
field measurements contains contaminants, dissolved solids, charged ions (such 
as sodium), and other impurities. The consequence of this assumption is that 
the diffusive terms within the liquid mass balance equation are zero 

v ij,l = v u = 

Assumption 7^5: The liquid phase is assumed to be incompressible. This assump- 
tion is valid under moderate pressures and allows us to remove the liquid phase 
material time derivative of density from the liquid mass balance equation 

D l p l 



Dt 



0. 



In isothermal conditions the density of the liquid phase can be assumed constant 
in space and time. In the presence of thermal gradients, on the other hand, we 
presume that the density of the liquid phase is a function only of temperature 
given by the empirical model 

p\T) = 10 3 (1 - 7.37 x 1(T 6 (T - 277.15) 2 + 3.79 x 1(T 8 (T - 277.15) 3 ) 

(5.39) 



measured in kg/m (and where [T]—K). See Figure 5.1(a) 
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Assumption 77=6: The gas phase is assumed to be an ideal binary mixture of water 
vapor and inert air. There are most certainly more than two species in most 
practical gas mixtures, but here we are concerned with with the diffusion, evap- 
oration, and condensation of water vapor within the gas mixture. The other 
species are assumed to be non-reactive and are therefore all grouped together 
into the air species. We choose the mixture to be ideal so that we can take 
advantage of the ideal gas law. This is valid since (a) the gas pressures under 
most experimental considerations are close to atmospheric, (b) under Richards' 
assumption [62J [65], the bulk gas pressure doesn't vary much under most ex- 
perimental considerations, and (c) the temperatures under consideration aren't 
far from standard room temperature. The use of an ideal gas mixture will 
break down under higher pressures, higher temperatures, and possibly under 
high variations in temperature. 

Assumption 7^7: The gas-phase chemical potentials and densities are only func- 
tions of the relative humidity and temperature 

fjpi = ^i fa T ) pfn = p9j fa T y (5 40) 

We make this assumption based on the fact that at the pore scale we can easily 
convert between the chemical potential, the density, and the relative humidity. 
Furthermore, this allows for us to tie the gas-phase mass balance equation to 
experimentally measurable quantities such as the relative humidity. 

Just as at the pore scale, we define the macroscale relative humidity, if, via the 
saturated vapor density, p sa t, and the density of the water vapor in the mixture: 

P 9v = PsatV, (5.41) 

where p sat = p sa t(T) can be expressed through the empirical equation 

lhat = «*P (31.37 -M14.79/T- 7.92 xlO-T) x ^ (( . 42) 

98 



(see Figure 5.1(b) ). 

The chemical potential of the water vapor is defined through the ideal gas law 



as 



fj,9v =f Jl Bv+R*'T\n(\<p) 



(5.43) 



where A = p sa t/p* is a function of temperature from (5.42 ) and p* is atmospheric 
pressure. 

The reason we are calling this an "assumption" is that, strictly speaking, these 
relationships hold for the pore-scale chemical potentials and pressures. We are 
dealing with averaged (upscaled) quantities so we make the assumption that 
these quantities follow the same functional forms. It is known that the upscaled 
pressure, density, and chemical potential are not the same as the pore-scale 
pressure, so in effect we are defining the upscaled relative humidity through 
these relationships. 
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Figure 5.1: Densities as functions of temperature 



Under assumptions 4 and 5 on the liquid phase we reflect now on the choice of 
the form of Darcy's law for the liquid phase. In the absence of species it may not be 
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reasonable to use the chemical potential form and instead revert to the pressure form. 



Recall from equation (4.77) that the Darcy flux for a fluid with one species is driven 



by gradients in Gibbs potential and temperature. Recall also that the coefficient of 
the temperature gradient is the macroscale entropy. To side step the necessity of 
modeling the liquid phase entropy and Gibbs potential directly we use the pressure 
form of the Darcy flux: equation ( |4.66[ ). Given one liquid species, a rigid solid phase, 
two gas species (see assumption #6), and the assumption that the gas densities are 
functions of temperature and relative humidity (see assumption #7), the Darcy flux 
for the liquid phase can be written as 



e'Jl' ■ (s l v l ' s ) = -e l Vp l - n l m Ve l - n l ^Ve g + e l p l g 



dip 5 



.dip 1 



dp 1 



dp 1 



+ <+<^'-VEEv,. 



]=v,a 



dip 1 



8p9j 



■■ -e l Vp l - e (tt'w - tt^)) VS + e l p l g 

( a M 9 s ,W\ <9p' 



dp 1 



dp 1 J dT 



*V£ 



3=v, a 



dip 1 

3p9i 



dp 



Hj 



dT 



■VT + 



dp 



,9j 



dip 



■Vip 



-e l Vp l - e (tt'« - 7r'w) VS + e l p l g - e l C l T VT - e l C l v V<p. (5.44) 



The functions C T and C are implicitly defined by equations (5.44) and may be 



functions of any variable (s) from the set of independent variables for the Helmholtz 
Potential. The coefficient of the saturation gradient can be rewritten as 



l l I dip 1 dip 1 

££p \W~d& 



2e l p 



! ,dip l 



as- c ' c - 



(5.45) 



This coefficient function measures the changes in liquid energy due to changes in 

saturation while holding density fixed. The notation chosen for these coefficients 
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is meant to be descriptive; the subscript indicates the associated gradient and the 
superscript indicates the phase. 



Dividing both sides of (5.44) by e gives the simplified pressure, saturation, tem- 



perature, and relative humidity formulation of the liquid Darcy flux 
g . ( £ l v l ' s ) = -Vp l + p l g - C l s VS - C l T VT - Cj,Vp. 



(5.46) 



The first two terms on the right-hand side are the classical Darcy terms, and the 
functions C l s , C l T , and C l are, as of yet, unknown. All of these new functions measure 
cross coupling effects due to the presence of other phases. Thought experiments used 



to make sense of these new terms will be presented in Section 5.4.1.1 after a deeper 
discussion of capillary pressure. 

Under assumptions #1 - #7, the heat and mass transport system can now be 
written as: 
dS 



at 



V ■ {K! ■ [W + C l s VS + C l T VT + C\,Vip - p l g] } 
M (p l - p 9 ») 



P' 



(p l -p 9 ») 



(5.47a) 



OS 



do 9v 

-V- {p 9v R gv -[Vfi 9v -g]} 

- V ■ {p 9 "K 9 ■ [p 9v Vfi 9v + p 9a Vp 9a + p 9 r] 9 VT - p 9 g] } 
= -M (p l - p 9 ") (p l -p 9v ) 
= pc p f - V • (JC • VT) + ph + Le\ 



(5.47b) 



+ 



p b 



P 



"i- 



'<-!■ 



d 



+ 2reS + T—(7r 9 -n l ) 



dT 



eS 



-e(l-S)H 9 p 9 + e(l-S) ^ 



j=v,a 



p 9 > - T 



dp 



tij 



+ 



+ 



dT 
■ (e l v l > s ) 

T 9 Vp 9 + J2 [/** Vp 9 '] + p 9 c 9 p VT 



p9i 



n idp l \„ m Tdpi, 

9 ^ dT) + e l dT 



3=v,a 
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L j=v,a ^ 



V [e 9 p 



9 n 9^ 



(e 9 v 9 ' s ) 



(5.47c) 
This system of equations originated from mass, momentum, and energy conservation 
and was supplemented with constitutive forms of the rates of mass, momentum, and 
energy transfer. We used the incompressibility of the liquid phase to arrive at the 
fourth line of the energy equation. In the gas phase, the change in pressure with 
temperature is given via the ideal gas law: 

f = (££) T, (5.48) 

V {?*\ + (™)B£ (5 , 49) 



dT \M9 J \M9 J dT 

where M 9 is the molar mass of the gas mixture and R is the universal gas constant. 
In the liquid phase, the change in pressure with temperature is the ratio of isobaric 
and isothermal compressibilities of liquid water 



dpi / i dv i\ J / i gyi\ a 



j 



dT \V l dT J I V V 1 dp 1 J P r 

Recall that p l = p l {T), p 9 - = p 9v (p,T), p 9 > = //*(</?, T), e l = eS, e 9 = e(l - S), 
and 7] a = r] a {T). Furthermore, e a v a,s is the Darcy flux associated with the a— phase 



(see equation (4.76)) and the latent heat, L, is an empirically based function of tem- 
perature. Therefore, assuming that the enthalpy, internal energy, and the linearization 
coefficients are known functions of these same variables, equations (5.47a) - (|5.47c) 



can be seen as a closed system of equations in saturation (S), relative humidity (ip), 
and temperature (T). It remains to find relationships for the linearization coefficients, 
the cross coupling Darcy terms, the gas-phase entropy, the enthalpy, and the chemical 
potentials. In the next subsections we discuss dimensional analysis, functional forms 
of the coefficients, and further simplifications for each equation one at a time. 

5.4.1 Saturation Equation 
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In the liquid phase, the linearization constant, K^ , is a function of the ease in 
which fluid flows through the medium. This is known as the hydraulic conductivity 
of the medium. The hydraulic conductivity is also known to be a function of the 
permeability of the medium. In saturated (rigid) media this is considered constant 
(or at least a tensor), but in unsaturated media they are typically taken as functions 
of saturation. In the present case, a careful inspection of the units indicate that 

k, k 

K l = = = ^, (5.50) 

— W P l 9 

where k is the permeability tensor of the medium, k is the hydraulic conductivity 

tensor, and \xi is the dynamic viscosity JSJES]- Notationally "/i a " (with a superscript) 

will denote chemical potential, and "/i Q " (with a subscript) will denote dynamic 

viscosity. 

The permeability, k, is typically separated into a saturated permeability, k , and 

a relative permeability, n ra . The relative permeability is assumed to be a function of 

saturation and depends on whether a is the wetting or non- wetting phase [62]. There 

are several functional forms of K ra , but one of the more commonly used is that of van 

Genuchten 1791, 



K r [ — K rw — \^e) 1 1 


"l - (S e ) 1/m 


m-i 2 


(5.51a) 


K r g i^rnw [J- WeJJ 


~l-(S e ) 1/m 


2m 

) 


(5.51b) 



where m is a fitting parameter, and S e is the effective saturation defined by 

£e = / ~_fr See [0,1]- (5.52) 

*-' max ^min 

Typical values of m are less than 1 where m = 2/3 is commonly used as a starting 

point for fitting numerical models to experimental data. Typical relative permeability 

curves are shown in Figure |5.2| The reader is to keep in mind that there are several 

such models in the literature (5j |62] . The van Genuchten model simply constitutes 

a widely used relative permeability model. Note that there is not a symmetry in 
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k rnw and k rw in the sense that k rnw (S e ) ^ k rw (l — S e ) as would naively be assumed. 
This is a manifestation of the fact that unsaturated media behave differently during 
imbibition and drainage. The value of k. is chosen based on the type of medium. If 
the medium is isotropic then the tensorial notation can be dropped and values from 
Table IE. 21 can be used. 

van Genuchten relative permeabilities (wetting=blue, nonwetting=red) 




o.:i 
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Effective Saturation 



0.8 



().!) 



Figure 5.2: van Genuchten relative permeability curves. The red curve shows the 
non- wetting phase, K rnw (S e ), and the blue curves show the wetting phase, K rw (S e ), 



each for m = 0.5, 0.67, 0.8, and 1. 



5.4.1.1 Capillary Pressure and Dynamic Capillary Pressure 

The capillary pressure, p c , is typically defined as the difference between the non- 
wetting (gas) and wetting (liquid) phase pressures when measured in a tube at equi- 
librium 



Pc Pnon— wetting Pwetting- 



(5.53) 
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At the microsale, the difference is related to the surface tension of the fluid, the contact 



angle, and the effective radius through the Young-Laplace equation (see Figure 5.3) 

27 cos 9 



Pc 



(5.54) 



The question is which pressure (thermodynamic, classical, or wetting (see Section 




Figure 5.3: Contact angle and effective radius in a capillary tube geometry. 9 is the 
contact angle, r is the effective radius, and k is the radius of curvature of the interface. 



4.3[ )) represents the non-wetting and wetting pressures in equation (5.53). The capil- 
lary pressure is measured with a force transducer in the same manner that the classical 
pressure is measured. For this reason we define the capillary pressure as 



Pc = P g ~P 1 - 



(5.55) 



Now that we understand which pressure is associated with the capillary pressure 
we turn to the entropy inequality to derive a constitutive equation equation for the 
time rate of change of saturation. In Richards' equation it is standard practice (as 



mentioned in Section 5.1.1 ) to take the capillary pressure as a function of saturation. 



These relations are reasonable for an equilibrium relationships. In the present model- 
ing effort we look toward the entropy inequality to determine an appropriate form of 



p c away from equilibrium. In the entropy inequality (equation (4.13)) there are two 
terms associated with the time rate of change of saturation: 



—pe and - p 9 e 9 . 
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Since e 9 = —e l these terms can be combined to give (p 9 — p)e . The time rate 
of change of volume fraction is a constitutive variable so the associated linearized 
equation is 



(p-f) ={p 9 -p 1 ) 



n.eq 



T€ 



eq 



(5.56) 



where the equilibrium state is not necessarily zero and the minus sign is chosen to be 



consistent with the entropy inequality. From the three pressures relationship, (4.52) 
the classical pressure is given as 

P = P + TV 



where p a is the thermodynamic pressure and ir a is a wetting potential. The difference 
in thermodynamic pressures is therefore rewritten as 

Pe ■= P 3 -P l = (p 9 - TT 9 ) - (p l - IT 1 ) =p c - (lT 9 - TT l ) =p c -TT c 



and equation (5.56) becomes 



(Pc ~ tt c ) 



n.eq 



(Pc ~ VT C ) 



TE 



(5.57) 



<</ 



Rewriting we get 



Pc 



Pc 



n.eq 



+ TT. 



<<l 



- 7T r 



n.eq 



eq 



— re 



(5.56 



We assume that the effect of the solid phase on the capillary pressure is completely 
captured by the preferential wetting, 7r c . Without the solid phase, the normal pres- 
sures of the liquid and gas phases are zero (this is the case with a flat interface). 
With this assumption the thermodynamic pressures are equal across the phases at 
equilibrium. Therefore, p c \ eq = 0. This implies that p c \ eq = p c \ eq + n c \ eq = 7i c \eq- 
Therefore the capillary pressure at equilibrium is interpreted as the difference in wet- 
ting potential and we arrive at an expression that is similar to that found in [47j. 
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To avoid possible confusion we will continue to use the symbols p c \ eq in place of n c \ eq 
even though they are understood to be the same. 

We finally arrive at an expression relating the classical liquid-phase pressure that 
appears in Darcy's law, p l \ n .eq, and the capillary pressure, p c \ e q- 



-P 



Pc 



n.eq 



+ 7T 



'-'-I 



T\< 



n.eq 



«/ 



\f 



re . 



n.eq 



(5.59) 



If the deviation in the wetting potential from equilibrium is assumed to be small 
relative to the pressure and the dynamic effects we can approximate the liquid pressure 

as 



~P 



Pc 



n.eq 



p y 



cy 



— re 



(5.60) 



n.eq 



where it is possible that p 9 ~ as well (in fact, this is a common assumption). To 
see why the deviation in wetting potential might be small, consider that in equation 
(5.58) if p c \n.eq ~ Pc\eq then the saturation dynamics is driven by the deviation in 



wetting potential. The deviation in wetting potential measures how much the shape 
of the curved liquid-gas interface is away from equilibrium. In slow flows it is unlikely 
that this deviation is significant. 



As mentioned in Section 5.1.1, the (equilibrium) capillary pressure can be related 
to the effective saturation through the van Genuchten p c — S relationship. This 
relationship depends on several fitting parameters and is given as 



Pc{S e 



g)(*^-i) 



1—771 



(5.61) 



where a has units of reciprocal pressure and m is the same fitting parameter as in the 



relative permeabilities (5.51a) j5j|62]. See Figure 5.4 for several examples of capillary 
pressure - saturation curves for various sets of parameters. Generally speaking, m 
increases (toward 1) as the soil becomes more densely packed. 
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van Genuchten p c — S curves 




o.i 



0.2 



0.3 



0.4 0.5 0.6 

effective saturation 



0.7 



0.8 0.9 



Figure 5.4: Examples of van Genuchten capillary pressure - saturation curves for 
various parameters. 

Substituting the capillary pressure and van Genuchten relationships into Darcy's 



law, (5.46), the fluid flux becomes 

R ■ (e l v l >*) = (^ " c£) VS - reVS + p'g 

- Vp° - CJ.V<f - C l T VT + V (vr c | n . eg - n c \ eq ) . (5.62) 

S v ' 

«o 
To understand the newly terms proposed here, we make the following three comments: 

1. First consider the gas pressure and relative humidity terms. In the absence 
of gravity, if the saturation, temperature, and the change in capillary wetting 



potential are held fixed then (5.62 ) states that flow is driven by gradients in rel- 



ative humidity and gas-phase pressure. The gas-phase pressure and the relative 
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humidity are proportional to each other where the constant of proportionality 
is a function of temperature and the species densities. With this in mind, these 
two terms together can be rewritten as a gradient in gas pressure. While a gra- 
dient in gas pressure can certainly cause flow, it is commonly assumed that p 9 
is approximately constant (known as Richards' assumption [62J) and therefore 
these terms are typically neglected. If these terms are not neglected then they 
are best written as a single gradient of relative humidity for easy coupling with 
the gas-phase diffusion equation 

2. Next consider the gradient of temperature term. In the absence of gravity, 



if saturation and relative humidity are held fixed then (5.62) states that flow 
is driven by a gradient in temperature. Saito et al. [67] indicated that the 
thermally induced flow was negligible as compared to isothermal flow (also 



discussed in [78, 80J). This indicates that the VT term in (5.62) is likely quite 
small. 

3. Finally we discuss the role of C l s = 2p l -^r. This function (or constant) relates 
the changes in energy with respect to saturation. The term is already associated 



with the gradient in saturation as seen in equation (5.62). From the VS" term 



in this equation we can see C l s as an enhancement of the capillary pressure - 
saturation relationship that directly models the affinity for the liquid phase to 
the other phases. It is entirely likely that this term is so closely linked with the 
capillary pressure that in experimental settings it is impossible to discern this 
effect from others. 

The saturation equation can finally be written as 



dS „ 

£ Tt~ v - 



K(S) ■ (l~ + C l s \ VS e + reVS e + C l T VT + Cj^ip - p l g\ 
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= *i (J - pf) tf - p») , (5.63) 

where we have assumed that 7i c \ n , eq — n c \ eq m and, abusing notation slightly, the C l 
term has be redefined to incorporate changes in the gas pressure. 

To account for the residual (minimum) saturation, the saturation is scaled to 
the effective saturation according to S e — (S — S min ) / \S max — S min ). Defining es 
as the product of porosity and the difference in maximal and minimal saturation, 
es '■= £(S max — S m in), an d letting S notationally stand for S e allows us to write the 
saturation equation as 

OS 



M l 



£sP l 



Ss'KiS) ■ {{-% + C s\ V£ + resVS + C l T VT + CjVp - p l g\ 
(p l -p 9v ) (^-/i 9 ") (5.64) 



At a quick glance, the sign of the VS term looks suspicious as it seems to indicate 
a backward heat equation. Observe that p' c (S) < for all values of S. Taking 
only the first line with C l s = returns Richards' equations exactly. The VS* term 
(henceforth referred to as the dynamic saturation term) was originally proposed by 
Hassanizadeh et al. in several publications (examples include [171111]) and is gaining 
more widespread acceptance in the porous media community. Taking all of the terms 
on the first line (again with C l s = 0) along with the dynamic saturation term gives 
a closed pseudo-parabolic equation in saturation. The C l s , C l T) C l terms along with 
the form of the right-hand side are all novel to this work. The temperature and 
relative humidity coupling terms can certainly be taken to be zero in certain physical 
instances, but generally the relative weight and functional forms of these terms is, as 
of yet, unknown. 

We now turn out attention to the gas phase diffusion equation. Analysis and 
numerical solutions to the saturation equation will be considered in Chapter [7j 

5.4.2 Gas Phase Diffusion Equation 
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In this subsection we make certain simplifications to the gas-phase diffusion equa- 
tion so as to tie the chemical potential formulation to the more classical enhanced 
diffusion model. As a first step toward this simplification we consider the fact that 



the gas phase chemical potentials are related to each other through equation (4.85); 
the expression for the relative motion of diffusing species in a binary system: 



N 

E 



F 



R9jT 



D a -[V^-g] 



0. 



With this, the gradient of chemical potential of the inert air in (5.47b) can be rewritten 
as a function of the water vapor chemical potential 



,9oT7„Sa 



p 9a Vp 



R9v 



(Vp 9v -g) + p 9a g. 



This means that the gas-phase mass balance equation can be rewritten as 
9 (ep£Ml - S)) 



dt 



V- lp 



-J-lv 



D 9v + p 



R9v 



K 9 



[S?(jL*>-g] 



v • {fffPrfK 9 ■ vr} = -M (p l - p 9v ) (fi l - p 9v ) . 



(5.65) 



Typically, one would choose a functional form of D 9v to match the enhance- 

and the functional form of K 9 from the van 



ment model discussed in Section 



5.1.2 



Genuchten model discussed in Section |5.4.1| In the present case we argue to use dif- 
ferent functional forms of D 9v and K 9 . This is done by considering the conversions 
between the pore-scale density and chemical potential to the relative humidity. For 
simplicity the tensorial notation is dropped and we assume that the diffusion and 
conductivity tensors are all scalar multiples of the identity matrix. 

We begin with some logical considerations for the gas-phase diffusion coefficient. 
If the gas-phase volume fraction were to drop to zero then there would be no gas in 
the pore space (or their would be no pore space) and the diffusion coefficient should 

drop to zero. Similarly, if the gas-phase volume fraction were to increase to 1 (100% 
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gas with no solid or liquid), then the diffusion coefficient should return to the Fickian 
diffusion coefficient D 9 . With these two limiting cases in mind we first propose that 
D 9v = Ce 9 D 9 where C is a scaling parameter. 

As seen in Chapter [2j the diffusion coefficient is modified for Fick's law based on 



the dependent variable of interest. In equations (2.1) and (2.3) we see a scalar factor 
of l/(R 9v T) between the mass and chemical potential forms of Fick's law. Making 
the same modification here along with the factor of e 9 suggested above we get 

p 9 *D 9 "V^ -► p* C^j D 9 Vp 9 " = ( £psat ^~ S) ^J D 9 V^ (5.66) 

where D 9 is the same pore-scale diffusion coefficient as found in Chapter [2} One 
simple way to look at this conversion is that it scales out the units and magnitude of 
the chemical potential when converting to relative humidity. That is, .D Su V p 9v and 
D 9 /(R 9v T)'Vip have the same units and magnitude. A further justification of this is 
found by recalling the pore-scale definition of the chemical potential: 

= p 9v +R 9v T In (\(p), (5.67) 

where A = p g s v at /P g an d ptat * s the partial pressure of the water vapor under saturated 



conditions. Taking the gradient of (5.67) and neglecting the temperature variation 

gives 

R 9v T 

VfJ, 9v « Vif. 

Hence we see the exact conversion used in Fick's law. 

Next we turn our attention to the hydraulic conductivity term that arose from 
Darcy's law: p 9v K 9 ^ p 9v . Similar to that of Fick's law, we need to scale the conduc- 
tivity to account for the fact that we're using the chemical potential as the dependent 
variable. Unlike the Fickian diffusion coefficient, this term already has the proper 

units since the units of p 9v ^ p 9v are the same as the gradient of pressure. Therefore 
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we seek a scaling that is unitless but scales the magnitude of the chemical potential 
down to that of pressure. That is, we need a constant, c, such that cp 9v K 9 'Vp 9v and 
K 9 "Vp 9 have approximately the same magnitude. 



Taking the gradient of both sides of the first line of equation (5.67) we arrive at 

'R9vTp 9 ' 



V> 



<-)v 



p 



3v 

, V ' " 
R9v Tp a\ f f i 



j)9v 
R 9v T 
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Vp 



9v 



Vp 9v 

R 9v T 



(.P 



Vp 9 . 



n9v 

— 1 Vp 9 



p9v I \ p9 

The coefficient of the gradient of gas-phase pressure can be rewritten as 

R°«T _ psatRQ-T _ p 9 / at _ A 

P 9 PsatP 9 PsatP 9 Psat 

Since the chemical potential form already has a factor of p 9v = p 9 s v at <p we scale K 9 by 
A to account for the difference in magnitude between the chemical potential and the 



pressure. Hence, the Darcy term in equation (5.65) is rewritten as 

R 9 °' 



P 



,!!•■ 



K 9 Vp 9v -> p 9v [ 1 



Rs* 



(XK 9 )Vp 9 \ 



Rgv j r y R9v 

Keep in mind that this is a scaling of the hydraulic conductivity; just as the factor of 
l/(R 9v T) is a scaling of the diffusion coefficient in Fick's law. 

One point of interest for this choice of scaling factor is that it is invisible when 
we consider a pure gas phase. That is, A = 1 when no species are considered since the 
saturated partial pressure will simply be the bulk pressure. This indicates that we 
have not actually changed Darcy's law. Instead we have simply made a conversion to 
account for the use of a different dependent variable. 



Next we focus on writing the gas-phase diffusion equation (5.65) in terms of 
relative humidity, saturation, and temperature. To do this we replace the chemical 



potential with relative humidity and temperature via equation (5.67). Taking the 



gradient of the chemical potential in equation (5.67) we get 



R 9v T 
Vp 9v = Vy2 



R 9 "T dX 
A dT 



+ R 9 " ln(A^) VT. 
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With the Fickian and Darcy terms written in terms of the relative humidity, along 
with the fact that the saturated vapor density is a function of temperature, the vapor 
diffusion equations can be written as 




at 



(ep sat (p(l - S)) 



V • <^ p sat (p 






if \ A dT 

- V • {p 9 p sat <PV g (Ai^ 9 ) VT} = -M (p l - p sat <p) (p l - n*) . 

Combining like terms, dividing by the porosity, replacing the hydraulic conductivity 
by the saturated and relative permeabilities, and simplifying gives 
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-V-{p sai 2%,S,T) 



Vip 
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V-{ Psat N9{<p,S,T)VT} 



M l (p l - pH") , , 



(5.68) 



where the functions T> and N 9 are 



2%, S, T) := (1 - S)D 9 + p 3at <pR°°T 1 



N°{<p,S,T):=<p 



R 9a 

Rgv 
1 d\ R* ln(Ay) 



ep g 



K rg (S) and (5.69a) 



+ p 9 Psatrj 9 ( ) Krg(S) 

ep g 



(5.69b) 
The enhancement model suggested by de Vries, and subsequently used by sev- 
eral authors [Ml E7J [75], EHJ ED], is a multiplicative combination of the pure Fickian 
diffusion coefficient, D 9 , the tortuosity, r = t(e 9 ), and an enhancement factor, r\: 



D = TT]D 9 . 



(5.70) 



In these works, the functional form of the enhancement factor is taken to be of the 
form suggested by Cass et al. 

e l 



V{a) 



a + 3- 



£ 



(a — ljexp < — 



1 + 



2.6 \ e 



fcj £ 



(5.71) 
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Here, f c is the mass fraction of clay in the soil. In the absence of clay the enhancement 
factor is taken as 



rj {a) = a + 3 



(5.72) 



(for an example where f c ^ see Saito et al. [67J). The tortuosity is taken to be a 
function of the volumetric gas content, 



t = (2/3)e 9 . 



(5.73) 



Using equations (5.72) and (5.73) in the multiplicative expansion of the diffusion 



coefficient, (5.70) gives a diffusion coefficient of 



D 



a + 3- 



D 6 



(5.74) 



The tortuosity and the porosity communicate to the diffusion coefficient the type of 



geometry under consideration. The present model (equation (5.68)) communicates 
this information via the porosity, the relative permeability, and the saturated perme- 



ability. The diffusion model using equation (5.74) relies on a fitting parameter, while 



the present model avoids this trouble. In the author's opinion, this highlights the 
main advantage to using the chemical potential as a modeling tool. 

Comparing the enhancement model of Cass et al. (using the material parameters 
from the experiment by Smits et al. [78]) to the present model, we see, in Figure 



5.5 that the relative humidity level curves of the present model underestimate the 
enhanced model for many values of the fitting parameter, a. That being said, these 
curves do suggest an enhancement over regular Fickan diffusion and, depending on the 
parameters of interst, give similar levels of enhancement as the model used in [75] . We 
simply state here that the present model offers a modified view of the enhancement 
model. There are several parameters that play roles in this model, but the advantage 
to the present approach is that all of the parameters are readily measured for a given 
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Comparison of Diffusion Coefficients 
Genuchten parameter: m=0.9438 (n=17.8), saturated permeability: k=1 .04e-f0 
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— Present Model 
— Enhancement Model 
--- Fickian Diffusion Model 
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Saturation 



Figure 5.5: Comparison of different diffusion models at constant temperature (T = 
295. 15K). The value for the saturated permeability was chosen to match that of 
[ks = 1-04 x 10 _10 m 2 ), where they found a fitting parameter a = 18.2. The 



"Present Model" refers to equation (5.68) (with VT = and no mass transfer) and 
the "Enhancement Model" refers to equation (5.3) along with (5.70), ( 5.72| ), and 
(5.73) for the diffusion coefficient, enhancement factor, and tortuosity respectively. 



medium (at least in laboratory experiments). There is no fitting parameter, so the 
type of material should dictate the level of enhancement. 

Another way to look at the present model is to consider that in most classical 
situations the gas-phase pressure is considered constant. The effect of this assumption 
is that the Darcy terms in the gas-phase mass balance equation are neglected. This 
assumption is valid in many cases, but in the present case the Darcy term is broken 
into component parts (air and water vapor) via the chemical potentials. The chemical 
potential formulation draws influence from the Darcy-type movement, along with 
the Fickian diffusion, of the individual constituents to define the general diffusion 
coefficient. 
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It is emphasized here that the traditional (de Vries-type) view of diffusion in 
porous media is not taken here. Shokri [73] suggested that the mechanism of enhanced 
diffusion is driven by the coupling of Darcy and Fickian diffusion. The novelty here 
is that the advection and diffusion are modeled in terms of the same dependent 
variable; the chemical potential. This suggests that the enhanced diffusion problem 
can be modeled by coupling Darcy-type flow along with Fickian diffusion in the gas 
phase. The relationship between the enhancement model and the present model will 
be discussed when we consider numerical solutions in Chapter [7| 

5.4.3 Total Energy Equation 

Continuing with the equation-by-equation derivation of the total heat and mois- 
ture transport model, we now turn out attention to the total energy equation. This 



picks up from equation (5.38) and we apply the simplifying assumptions presented in 



the beginning of Section 5.4 

If we assume that the vapor and air densities are functions of relative humidity 



and temperature only, the total energy equation (5.38) can be written as 
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(5.75) 
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Recall that p = p{<p,S,T), Pc = p c {S e ), p?> = n*{<p,T), p» = p*{<p,T), T 9 = 
T 9 (ip,T), rf = r] 9 (T), p l = p l {T). Also recall that e a v a,s represents the Darcy flux 
for the a phase: 



£ l v l,s = _ K l 






-K 9 



{-p' c {S e ) + C l s ) VS e + reVS e + C l T VT + C\,V<p - p l g (5.76a) 
R 9a 



Xp 



,!-)>• 



1- 



R9v 



d ~W VT + W Vlf ' + pVVT " p ^ 



(5.76b) 



It is clear that there are several physical processes and couplings that occur for 



energy balance to be achieved. Equation (5.77) below shows the classical 1958 model 



of de Vries [32] (which is similar to that of Bear [3] and is also presented in [T3] for 
the saturated case). 



" C 'f - * (P'W - >W) f 



V-(|£Vr)-Le<-(Wfn 



f e « v «,^ I VT. 



(5.77) 



In this form of the energy equation, W a is a differential heat of wetting [2], and 
the other variables are written in the present notation for convenience. At first 



observation, the T, S, K^, e l , and VT terms in equation (5.75) are similar to terms 



found in the de Vries model. That is, we capture the standard effects of specific 
heat along with differential heat of wetting, thermal conductivity, mass transfer, and 



convective heating. Implicit in the S term in (5.75) is that we relate the partial 



derivative of the difference in wetting potentials, Td(ir 9 — ir l )/dT, as a differential 
heat of wetting. The present model also captures the effects of changing relative 
humidity, nonlinear effects such as VS 1 ■ V5 and V<y2 ■ V<£>, and cross effects such 
as VS 1 ■ V(/?. It remains to determine which (if any) of these effects are negligible 
as compared to the others. To make this determination we perform a dimensional 
analysis in the next subsection. Let us first focus on the thermal conductivity term, 

V • (K ■ VT) . 
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The functional form of the thermal conductivity, -ff , can be approximated in 
several ways. A first approximation is to take the thermal conductivity as a weighted 
sum of the conductivities of the individual phases 

a 

Comparing to results in [77], we note that this seems to overestimate the measured 
thermal conductivity as well as fail to capture the experimentally measured curvature 
of the thermal conductivity - saturation relationship. Since K^ is a linearization 
constant that arose from the entropy inequality, it can depend on any variable which 
is nonzero at equilibrium. In particular, Ki is a function of saturation. Smits et 
al. [77] use a combination of the Cote-Konrad and Johansen models to estimate the 
thermal conductivity in the scalar case: 

K(S) = K e (S) (K sat - K dry ) + K dry , (5.79) 

where K sat is the conductivity of the saturated medium, K dry is the conductivity 
of the dry medium, and K e (S) is a "normalized thermal conductivity known as the 
Kersten number." Cote and Konrad proposed a functional form of K e as 

*.W = Y^^s- (5-80) 

The parameter, k, is a fitting parameter that is presumed to be different for each type 
of soil. In [77] , k was estimated for several types of sands and several types of soil 



packs. Figure 5.6 shows a thermal conductivity curve for (5.79) with tightly packed 
30/40 sand that has a porosity of 0.334. For comparison, equation (5.78) is shown in 
red for the same experiment. 

We make some comments here giving some possible reasons for the discrepancy 



between the weighted sum model (equation (5.78)) and the model that more closely 



matches what is experimentally observed (equation (5.79)). First, the thermal con- 



ductivity of air is neglected as compared to the thermal conductivity of water or 
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Comparison on Johansen-Cotc-Konrad and Weighted Sum Thermal Conductivitcs 
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Figure 5.6: Johansen thermal conductivity model with Cote-Konrad K e — S relation- 
ship (with k = 15) plotted in blue, and the weighted sum of the thermal conductivities 
of the individual phases plotted in red. 



solid. Also, the thermal conductivity of liquid is much smaller than that of the solid, 
K l < K s . Furthermore, the geometry of the packed solid plays a crucial role. Observe 
that if we idealize the soil grains as individual spheres then there are relatively few 
contact points between the individual grains of the solid phase. This idealization can 
be used as a partial explanation for the left-hand tail seen in the Cote-Konrad model 



depicted in Figure |5.6| If there are few contact points between the individual grains 
then it is much harder for heat to transfer in the absence of a liquid phase connecting 



them. Thus, equation (5.79) tells us that all pertinent information is obtained by 



knowing what the thermal conductivity of the dry and saturated porous media is as 
well as an interpolation function for effective saturation. This captures more of the 
microscale geometry than just the volume fractions. The effects of these two proposed 
thermal conductivity functions on the behavior of the heat transport model will be 



explored when we consider numerical solutions in Section 7.4.2 
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5.4.3.1 Dimensional Analysis 

To determine which, if any, terms can be neglected from the energy transport 

equation we perform a dimensional analysis. Begin by noting that K/(pc p ) has units 
of area per time. This suggests a natural choice of time scale for the thermal problem 
of 

where t' is dimensionless time. Dividing by pc p (measured at a reference state), 
introducing x c as a characteristic length (e.g. the height of a column experiment), and 
multiplying by t c = (pc p x 2 c )/K gives the dimensionless form of the energy equation 
(the statement of which is suppressed for the sake of brevity). 

Recall that the volumetric heat capacity, pc p , is linearly related to the specific 
heats of the individual phases 



pCp = J2 (e a P a c a p ) = £S P l i + < l - S )P 9c P + (1 - £ )P 1 



v 



Taking S = 1 as a reference state (or equivalently, S = 0) gives a characteristic value 
of pc p . Using values from Appendix [El we see that pc p ~ O(10 6 ). Hence, several of 



the quantities in (]5.75|) can be neglected: 
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(5.81e) 
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In order to make these approximations it is assumed that Gibbs potentials are given 



by the Gibbs-Duhem relationship, (4.61), and that the Helmholtz potential and in- 



ternal energy are approximately the same order of magnitude as the Gibbs potential. 
With these considerations we can rewrite the present version of the energy equa- 
tion as 
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(5.82) 



(1 - S) dT 

Unfortunately this analysis leads us to the conclusion that this new version of the 
heat transport equation is only slightly different than those proposed in past works 
[Hl|32]. The major differences are the VS" terms associated with the Darcy fluxes, the 
capillary pressure adjustment to the differential heat of wetting term, and the Darcy 



fluxes themselves. Recalling the forms of the Darcy fluxes from equations (5.76), the 
energy equation can be rewritten in a more compact notation as 



dT 
=pc p — - V • (K ■ VT) + ph + Le l g + W 

+ {xiVS + X2 VT + X3 V^) • VT 

+ (x 4 VS + XsVip) ■ Vip + x 6 VS • VS 



as 

~dt 



(5.83) 



where W and each Xj are implicitly defined via equations (5.82) and (5.76). It remains 



to determine the functional form(s) of the several constitutive variables in (5.83). 



5.4.4 Constitutive Equations 



Hidden within the coefficients of (5.83), (5.64), and (5.65) are a few final re- 



lationships necessary for closure. In particular, we need constitutive equations for 



dpc 
de l 



(5.84a) 
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ef =M(p l - p 9 *) (/J - //") (5.84b) 

W = Pc(S) + 2tsS + T— {ix 9 - u l ) = Pc (S) + 2tsS + W (5.84c) 

C l s = 2p 1 ^ = e (ir l w - 7T 1 ^) (5.84d) 

The simplest possible assumption would be that r, C l s , C l T) C l , and W are constants. 
This would allow for the easiest sensitivity analysis but is likely contrary to physical 
reality. The following paragraphs discuss each of these terms and propose functional 
forms in terms of saturation, relative humidity, and temperature. The sensitivity of 
the numerical solution to several of these parameters is discusses in Chapter [7j 



It is generally assumed that r in equation (5.84a) is constant [HI EO], but ac- 
cording to the linearization process in HMT, r can be a function of any variable that 
is not zero at equilibrium. In particular, it is possible that r is a function of S; 
but which function? In [T7], the authors suggest several functional forms (constant, 
linear, quadratic, Gaussian, and error) and compare to experimental findings. Their 
findings suggest that ". . . an error function or Gaussian relationship for the damping 
coefficient r provides reasonable agreement between data and simulations." Thus we 
consider the following forms: 

r = T max (5.85a) 

r = W A _ erf (tLzJi\ \ (5.85b) 

r = r max exp I — — J . (5.85c) 



Plots of equations (5.85) are shown in Figure 5.7 with typical mean and standard 



deviation parameters. To the author's knowledge, no other experiments have been 

conducted to make a better determination as to the functional form of r. This 
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being said, since r is a measure of the rate at which the pore-scale saturation profile 
rearranges in a dynamic situation, it is reasonable to assume that as S — > 1 the 
effect of this term should be minimized and as S — > the effect should be maximized. 
Hence, in the author's opinion an error function is more sensible. It remains, of course, 
to determine the values of the maximum, mean and standard deviation parameters 
which are likely themselves functions of material properties. 

Proposed functional forms of r 














— error function 

— Gaussian 

— - - constant 






/ \ \ 




- 



0.1 0.2 0.3 0.4 0.5 0.6 0.1 

effective saturation f-1 



0.8 0.9 



Figure 5.7: Three proposed functional forms of r = t(S) 



The evaporation rate term, ef", given in equation (5.84b) is written as a function 



of the difference between the liquid and vapor chemical potentials. The chemical 
potential in the water vapor is a function of temperature and relative humidity 



fj,Bv =ij*>+R*>Thi.(\<p). 



The liquid-phase chemical potential, on the other hand, does not have such a natural 
description. At equilibrium, p) = fi 9v . Away from equilibrium we only know that 

y 
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/' 



r l = ^ l + p - 1: 



and therefore is a function of every variable that ip l depends. In the most simplistic 
form we can assume that the liquid chemical potential is yr = p[ + (p l — p{)/p l - This 
assumption is taken from classical thermodynamics (see [21J for example). Further- 
more, /4 « /if" if we take the reference state to be equilibrium. Therefore, 

ef « ^ ( P l - f) (i^. - R»Tln (A*,)) 

- ^ ( P l ~ p») ( ~ PC + T f~ pl ° ~ &-T In {\<p)\ , (5.86) 

where M is a fitting parameter. The factor of relative humidity is included to achieve 
a better match with existing empirical models (discussed in the next paragraph). 

There are several empirical rules for evaporation in porous media. One such rule, 
given by Bixler [Tj5] and repeated in Smits et al. [75], is 

ef = b(e l - e l r )R^T (p sat - p*) , (5.87) 

where b is a fitting parameter and e l r is the residual volumetric water content. Equa- 



tions (5.86) and (5.87) are quite different, but under proper scaling they are close as 



seen in Figures |5.8| From these plots it is also clear that there is a large dicrepancy 
between these model at very low saturations. These plots are generated at stan- 
dard temperature with S = 0. The dynamic saturation term will change the shape 



of these curves, but as the Bixler model, (5.87), is not dynamic we compare only 



with the steady state form of (5.86). Furthermore, the present model depends on 



the van Genuchten parameters for capillary pressure. In Figures 5.8 the parameters 
m = 0.944 and a = 5.7 are used along with b « 2.1 x 10 -5 to match the values used 
in [78]. 



The differential heat of wetting, W, in equation (5.84c) represents the heat gained 
or lost due to changes in saturation and adsorption. The present generalization 
suggests that the differential heat of wetting be supplemented by the capillary pressure 



and time rate of change of saturation. According to [64], the typical value of the 
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Level Curves of Evaporation Rate 



Level Curves of Evaporation Rate 




0.1 0-2 0-3 0.4 0.5 0.0 0-7 0-8 0.9 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Effective Saturation Relative Humidity 

(a) Comparison of mass transfer rate vs. ef- (b) Comparison of mass transfer rate vs. rela- 
fective saturation shown with level curves in tive humidity shown with level curves in satu- 
relative humidity. ration. 

Figure 5.8: Level curves of mass transfer rate functions. 

differential heat of wetting is on the order to 10 3 J/ kg depending on the type of soil. 
This value will be taken as constant throughout, but in reality value should be a 
function of saturation. 

Finally, the values of C l s ,Cip, and C l in equations (5.84d) - (5.84f) are new and 



hence there is no existing literature for which to make estimates or comparisons. 
For this reason we make the initial assumption that these terms are constant. This 
allows for relatively simple sensitivity analysis without introducing any unnecessary 



mathematical difficulties. As discussed in Section 5.4.1.1, the value of C l T is likely 
quite small since some research has been done to determine the affect of thermal 
gradients on Darcy flow |67j . 

5.5 Conclusion and Summary 

In this chapter we have derived several new equations and terms for heat and 
moisture transport in unsaturated porous media. For the sake of readability, we 
summarize the results, assumptions, and equations derived here within Chapter [5] 

The main assumptions are: 



Assumption #1 The solid phase is rigid, incompressible, and inert. 
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Assumption 77=2 The liquid and gas phases are composed of N constituents, (this 
was later relaxed to let N = 2 in the gas phase and N = 1 in the liquid phase). 

Assumption 77=3 No chemical reactions take place in any of these phase. 

Assumption 77^4 Diffusion with the liquid phase is negligible compared to the ad- 
vection of the liquid phase. 

Assumption 77=5 The liquid phase is incompressible. 

Assumption 77=6 The gas phase is an ideal binary gas mixture of water vapor and 
inert air. 

Assumption 77^7 The gas-phase chemical potentials and densities are functions of 
relative humidity and temperature. 

The secondary assumptions used up to this point are (in order of appearance): 

• the medium of interest is granular so angular momentum conservation yields a 
symmetric stress tensor, 

• the material is simple in the sense of Coleman and Noll [27] , 

• the phase interfaces are assumed to contain no mass, momentum, or energy, 

• second-order effects in velocity are negligible (e.g. v aj,a ® v aj,a <C v a ), 

• the species in the solid phase do not diffuse, 

• inertial terms in the momentum balance equation are negligible, 

• the capillary pressure - saturation relationship is given by the van Genuchten 
function, 

• the deviation in wetting potential is approximately zero ((7T c \ n . eq — n c \ eq ) ca 0), 
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the coefficient of the dynamic saturation term, r, is constant, 



Considering assumptions #1 - #7 along with all of the secondary assumptions, 
the final system of equations proposed to model heat and moisture transport in un- 
saturated porous media is: 

dS _ 



at 



e- s l K l {p' c + C l s ) VS + Ts s VS + C l T VT + Cj,V^ - p l g 



:-<-)u 



0_ 
di 



{p sat ip{l -S))-V- [ Psat {VVip + N°VT)] = -e 



<-lu 



dT dS 

= P c p — + W— -V ■ (K-VT) + ph + Le l g 

+ {xiVS + x 2 VT + XaVy.) • VT 

+ (x 4 vs + x 5 v^) • Vip + xe vs • vs, 



(5.88a) 
(5.88b) 



(5.88c) 



where the relevant empirical, constitutive, and derived relations are 



ft'.s 



K l (S) = -K rl = - 
Hi Hi 



H.i m ( 



yfs{\- [l-S 1 ^]^' 



K 9 {S) 



Hg 

1 



K rg 



1l (i _ S ) 1/3 (1 - s 1 ^) 

Hg 



2rn 



Pc (S) = - {s- 1 ^ - 1) 
dp c 



l—m 



t = — -r (see equations (5.85)) 
oe L ^ r 



f> = Mcp (p l -(?>)[ Pc + T f p o _ R^T In (\<p) 
V{<p, S, T) := i I - S)D* + Psat <pR 9 *T (l -?£)(^L] Krg (S) 



N9(v,S,T) = v 



R9 v J \ep g 
1 d\ R*> ln(A^) \ 
XdT + T J 



Sfig 



(5.89a) 
(5.89b) 

(5.89c) 
(5.89d) 

(5.89e) 
(5.89f) 



Vfa S,T)[{^=+ - Z X "^ ) + P'PsatV 9 ( — ) Krg(S) 



D 9 {T) = 2.12 x 10" 
e p a c\ 



T 



273.15 



(5.89g) 

(5.89h) 

(5.89i) 
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ph = J2p a h C 



(5.89J) 



W = p c + 2reS + iy 



, l + («-l)£ 

^(T) = 3.79 x 1(T 5 (T - 277.15) 3 - 7.37 x 10" 3 (T - 277.15) 2 + 10 3 

1 / 6014 79 \ 
p sa t{T) = -exp (31.37 - 7.92 x 1(T 3 T — J lO" 3 

fn(T) = (-2.56109 x 10~ 6 (T - 273.15) 3 + 0.00057672(T - 273. 15) 2 



Vg( T ) 



-0.0469527(T - 273.15) + 1.75202) 10" 3 
1.02312T 3 3.62788T 2 



— + 0.00665915T + 0.11767 ] 10" 5 

10 6 J 



10 9 
L{T) = 2.501 x 10 6 - 2369.2(T - 273.15) 

r] 9 (T) = 6.1771 x 10" 4 (T - 273.15) 4 - 7.3971 x 10~ 2 (T - 273.15) 3 

+ 3.1324(T - 273. 15) 2 - 34.4817(T - 273.15) + 191.208. 



(5.89k) 


(5.891) 


(5.89m) 


(5.89n) 



(5.89o) 
(5.89p) 

(5.89q) 
(5.89r) 



Equations (5.88) coupled with equations (5.89) give several adjustments to the 



classical models for saturation (Richards'), vapor diffusion (Phillip and de Vries), and 



heat transport (de Vries) presented in Section 5.1 In order for the present models 
to be accepted in the hydrology community we must show that the proposed terms 
are non-negligible and in some way put some of the empirical relations on a firmer 



theoretical footing. The proposed vapor diffusion equation (5.88b) is a prime example 



of this as there are no empirical fitting parameters within the diffusion coefficient 
(hence removing the need for an empirical enhancement factor) . 

In Chapter|6]we discuss the mathematical questions of existence and uniqueness of 
solutions to the individual equations. In Chapter [7] we discuss numerical simulations 
of the models. 
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6. Existence and Uniqueness Results 

In this chapter we discuss the necessary regularity and assumptions for existence 
and uniqueness of solutions for the three equations. As the main thrust of this work 
is not to prove existence and uniqueness for general classes of systems of partial dif- 
ferential equations, we approach these problems by stating relevant existing theorems 
from the literature and satisfying the hypotheses of these theorems. The saturation 
and gas diffusion equations are both of parabolic type and can be treated similarly. 
The heat transport equation is an advection-reaction-diffusion equation that, in prin- 
ciple, should be parabolic in nature. The advection terms force a different approach 



to this equation. In Section 6.1 an existence and uniqueness result for the saturation 



equation with the third-order term (due to Mikelic [57]) is outlined. The theorems of 



Alt and Luckhaus [U E] are outlined in Section 6.2 and then used in Sections 6.2.1 



and 6.2.2 to prove existence and uniqueness results for Richards' equation and the 



vapor diffusion equation respectively. Finally, an existence and uniqueness result for 



a special case of the heat transport equation is presented in Section 6.3 



6.1 Saturation Equation with r^O 

The saturation equation has been well studied since Richards' first introduced it 
in the 1930's. Recent modeling efforts, including those of Hassanizadeh et al., have 
introduced a new term into the classical Richards' equation and this has caused a 
resurgence in the analytical study of the saturation equation. The 2010 paper by 
Andro Mikelic [S7J gives the necessary conditions for existence and uniqueness of a 
weak solution to the following equation: 

f = V- {if (S) (-§V5 + rv(f )+«,)} in tt-nx (0,70 (6.1a) 

S = S D onT D = d D nx (0,T) (6.1b) 
K(S) (-^ VS + rV (^\ +e 3 Yu = RonT N = d N tt x (0, T) (6.1c) 

S(x,t = 0) = S l (x)onn (6. Id) 
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Here, e 3 is a unit vector pointing the z direction to account for gravitational effects, 
v is an outward pointing normal, and the subscripts D and N represent Dirichlet and 



Neumann conditions repsectively. Notice that (6.1) is a simplification of the present 
saturation equation as it contains no evarporation (source) term and no coupling with 
relative humidity or temperature. 

Mikelic's theorem is stated here for completeness. 

Theorem 6.1 (Mikelic 2010 |57| . Theorems 3 & 4) Consider the following hy- 
potheses: 

HI: there are constants /3 > 0, Sk > and a nonnegative function f G Cq°(R) such 
that K is given by 

*<-•> = i + a^/i-) -' 6 ' - 1 ' 

H2: there exists A > 0, S p > 0,M p > and an arbitrary function g G C£°(M) such 
that —p' c is written as 

H3: the product of the functions K and p' c is bounded on [0, 1]. 

H4: the initial Dirichlet data is smooth: Sp G C 1 ([0,T]; H l (Q)), and is bounded 
away from zero < SD,min < So(x,t) < 1 (or impose that Sd ^ Oa.e.) 

H5: R = i?oC; where R G C 1 ^ x [0,T]), R > and ( G Cg°(R), C(z) > /or z > 
0, and z£(z) > /or z < 0. 

H6: Initial moisture content satisfies a "finite entropy" condition: J Q (So(x)) ~ dx < 
+oo where (3 > A > 2 



Under these hypotheses there is a weak solution for (6.1) where S G H 1 ^?) such that 
0<S(x,t)a.e. onVt T , Vd t S G L 2 (fi T ) and 5 - S D GL 2 (0,T;1/) /or F = H l (0, 1). 
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The proof of this theorem is beyond the scope of this work, but it indicates 
that under constant relative humidity and temperature conditions, where no mass 
transfer is expected, there exists a weak solutution to the saturation equation. The 
sixth hypothesis restricts the shape of the initial condition. Simply put, the initial 
condition cannot drop to zero in such a way as to make J n S ~ dx go to infinity. This 
avoids the natural degenerate nature of the problem. The regularity expected for 
the solution (H l ) is a nice result given that this is actually a third-order differential 
equation. 

6.2 Alt and Luckhaus Existence and Uniqueness Theorems 

We now turn our attention to demonstrating the necessary conditions for exis- 
tence and uniqueness of Richards' equation (saturation with r = 0) and the vapor 
diffusion equation in the special cases where the other dependent variables are held 
fixed (possibly even constant). The two equations are treated together in this section 
since they both fall under the class of quasi-linear parabolic equations. As such, they 
can be analyzed using similar theory. For the purposes of demonstrating existence and 
uniqueness we apply general theorems by Alt and Luckhaus [1, 2 J to these equations. 

The following paragraphs are paraphrased from Alt and Luckhaus [2] and are 
presented here to introduce the reader to the notation used therein and for future 
reference. 

Consider the general initial boundary value problem (IBVP) for a system of quasi- 
linear elliptic-parabolic differential equations 

d t V{u) - V ■ o i (6(u), V«) = f j (b(u)) in (0,T) xfij' = l:m (6.2a) 

b(u) = b° on {0} x Q (6.2b) 

u = u D on(0,T)xr (6.2c) 

a j (b(u),Vu)-u = on (0,T) x (dn\T), j = 1 :m (6.2d) 
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In equations g2), u G R m , b : R m -> M m , and a : M m x M mxiV -> M m where JV is 
the spatial dimension of the problem and m is the number of equations. 

We call u in the affine space u D + L r (0,T; V) a weak solution of (6.2) if the 



following two properties are fulfilled: 

1. b(u) G L°°(0,T; L\Q)) and d t b(u) G Z/*(0,T; V*) with initial values b°, that is 



(d t b(u),C)+ [ f (&(«) - b°) d t ( = 

Jo Jn 

for every test function ( G L r (0,T; V) n W 1 ' 1 ^,! 1 ; L°°(fi)) with C(T) = 0. 



2. a(6(it), V«), f(b(u)) G Z/*((0, T) x Q) and tt satisfies the differential equation, 
that is, 

/ (d t b(u)X)+ f /a(6(u),Vu).VC= / [ f(b(u))( 

Jo Jo Jvt Jo Jn 

for every C G Z/(0,T; !/)• 

Recall from Functional Analysis that V* is the dual space of the vector space V, and 
W k, P = | w e L p^^ . ^^ e L p (n)\/\a\ < k} with the weak derivative D a u. The 
reader should also recall the common simplified notation W k,2 (Q) = H k (Q). 
Consider the following hypotheses: 

HI: Q C W 1 is open, bounded, and connected with Lipschitz boundary, T C dQ is 
measurable with if n_1 (r) > and < T < oo. 

H2: b is a monotone vector field and a continuous gradient, that is, there is a convex 
C 1 function $ : R m ->■ R with 6 = V$. We can assume that 6(0) = 0. The 
convexity of $ then implies that we can define 

B(z) := b(z) ■ z - $0) + $(0). 

H3: a(b(z),p) is continuous in z and p and elliptic in the sense that 

(a(b(z),p^)-a(b(z),pW)) ■ (p« _ p (2)) > c| p (i) _ p (2)| r 

with 1 < r < oo and f(b(z)) continuous in z. 
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H4: The following growth condition is satisfied: 

|o(6(*),p)| + \f{b{z))\ < c (1 + B{zf-^ r + Ipr 1 ) . 

H5: We assume that u D is in L r {0, T; W 1>r (fi)) and in L°° ((0, T) x Q) and we define 

V :={ve W^iQ) : v = on T} 

H6: Assume either that b° maps into the range of b and therefore there is a measur- 
able function u° with b° = b(u°) or that 

a t u D eL\0,T]L° o (£l)). 

The existence and uniqueness theorems of Alt and Luckhaus [2] are stated here 
for convenience and reference. 

Theorem 6.2 (Alt and Luckhaus |2], Theorem 1.7) Suppose the data satisfy 
HI - H6, and assume that dtU D G L 1 (0,T; L°°(fl)). Then there is a weak solution to 



(6.2). 



Theorem 6.3 (Alt and Luckhaus |2], Theorem 2.4) Suppose that the data sat- 
isfy HI - H6 with r = 2 and 

a(t,x,b(z),p) = A(t,x)p + e(b(z)) 

where A(t, x) is a symmetric matrix and measurable in t and x such that for a > 

A — al and A + adtA 

are positive definite. Moreover assume that 

\e(b(z 2 )) - e(b( Zl ))\ 2 + \f(b(z 2 )) - f(b( Zl ))\ 2 < C (b(z 2 ) - b( Zl )) (z 2 - Zl ) . 

Then there is at most one weak solution. 
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6.2.1 Existence and Uniqueness for Richards' Equation 

The existence and uniqueness of weak solutions to Richards equation is known, 
and a general tool for handling this problem is the Alt-Luckhaus theorem stated 
above. In this subsection we set up and state the theorems. We will show that the 



hypotheses of Theorems 6.2 and 6.3 are satisfied under restrictions on the Dirichlet 
boundary conditions and appropriate boundedness assumptions. The equation 

OS d fn rw (S)dp l 

K rw [o ) 



dt dx \ p l g dx 

= fa Ur W (S)— - K rw {S)j (6.3) 

is Richards' equation in dimensionless time and one spatial dimension. In this formu- 
lation we take h as the hydraulic head: h = p l /{p l g)- Recall from previous discussions 
that the pressure (or head) is a function of saturation. This relationship is invertible 
so here we note that saturation can be written as a function of pressure (or head). As 
suggested in [2JE3], "saturation may be less regular than pressure, therefore we expect 
to achieve better [regularity] results by applying a Kirchhoff transform" . A Kirchhoff 
transformation gives a smoothed relationship between head and a new unknown; a 
generalized pressure head 

JT : R ->■ R as X(h) = / K rw (S(q))dq := u. 

J — oo 

The pressure head is taken to be negative by convention (opposite sign of capillary 



pressure). The variable u now becomes the primary unknown of (6.3 ) since the spatial 
derivative can be written as 

du dJ(f dh dh 



dx dh dx rw dx' 
and if b(u) is defined as b(u) = S(J%^~ l (u)) — 1 then 



!«")) = |(|— ««)+!))■ (6-4) 
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Notice that the definition of b depends on the invertibility of JfT. Also, since 
<J>ff~ l (u) = h, b can be seen as a function of h: b(h) = S(h) — 1 (see Fig- 
ure 



6.1) . Given the van Genuchten capillary pressure - saturation relationship, 



S(h) = ( (ah) + 1 ) , and the van Genuchten relative permeability function, 

K-rw(S) = y/S (l — (l — S 11 /" 1 ) ) , Figure 6.2 shows several plots of J^ for various 



parameter values. There is a horizontal asymptote as h — > — oo and it is evident from 
the plot that J^ is one-to-one and onto for all values of h 6 (— oo, 0), but as h gets 
large in absolute value the inverse becomes unstable. 

b(h) = S(h) - 1 for m = 0.8 




Figure 6.1: The function b(h) = S(h) — 1 for m = 0.8 and various values of a. 
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Kirclilioff transformation for m = 0.< 



I). ir, 




5 - 10- 2 - 



Figure 6.2: Kirchhoff transformation J^ for m = 0.8 and various values of a. 



Matching to equation (6.2a) we note that j = 1, define a(b(u), Vu) as 



/ du \ du 

a f b(u), — )=— + K rw (b(u) + 1), 

and notice that / = 0. Given constant head Dirichlet boundary conditions, we finally 
rewrite Richards' equation as 



I du \ 

Wt {h{u)) = Tx\d~x~ Krw{h{u)) ) in (°' T ) xfi 

u = u D on (0,T)x{0,l} 
u = Uq on OxO 



(6.5a) 

(6.5b) 
(6.5c) 



With this form of Richards' equation we propose the following existence and unique- 
ness result. 



Theorem 6.4 Suppose that the following conditions hold for the generalized head, u, 



in equation (6.5). 
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1. tt = (0, 1) and t E (0, T) C (0, oo) 
£ m d GL°°((0,T) x tt) 
3. d t u D eL^O.T;! 00 ^)) 



T/ien i/iere ezzste a unique weak (u £ ud + L 2 (0,T; Hq(Q))) solution to (6.5) 



The proof of Theorem 6.4 has been discussed in several articles. In particular, the 



transformation of Richards equation to the form seen in equations (6.5) are discussed 
as a model problem for the Alt and Luckhaus theorems [2]. Furthermore, this proof 
is presented in [63] as part of their numerical formulation of Richards' equation. The 
fundamental reason for presenting this result here is that the value of r in the new 
saturation equations is not yet well known in experimental studies. Presenting this 
case simply covers all of the possible bases. 

In the cases where V</>, VT, or mass transfer terms are non-zero, these terms 
become source terms that depend on x. This means that / = f(x,b(u)) ^ 0. Ac- 
cording to section 1.10 in [2] "it makes no difference if a and / depend on x. n This is 
made more clear in their subsequent work [T] where the theorem is explicitly stated 
to allow for x and t dependence. 

For comparison sake we observe the difference between the regularity required for 



Richards equation (Theorem 6.4) and for the extended saturation equation with the 
third-order term (Theorem 6.1). For the equation with the third-order term (V-VS 1 ), 
the weak solution is in H l while in the second-order equation the weak solution is in 
L 2 . This extra required regularity is expected. 

6.2.2 Vapor Diffusion Equation 

To prove existence for the gas diffusion problem we proceed using the theorem 

of Alt and Luckhaus as with the saturation equation. Recall that under constant 
temperature and fixed saturation conditions 

(i-aw)^-s(^ s W)|) = ° (6 - 6) 
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where 

/ R 9a \ 
2%, S(x)) = (1 - S(x))D° + Psat! pR g »T 1 1 - — J AK(S(x)). (6.7) 

Allowing 5 to be a function of x constitutes a departure from the exact form of 



the parabolic-elliptic system found in Alt and Luckhaus (see equations (6.2)) as this 
is now a non-autonomous differential equation. In [lj this proof was generalized 
to allow for b = b(x,u) and for a = a(i,M,Vu) (see section 11 of PQ). The only 
additional assumptions for the existence theorems are that b : Q x R — y R and 
a : Q x R x R — >• R^ are measurable in the first argument and continuous in the 



others. With this addition to Theorem 6.2 we proceed with the existence theorem for 
the gas-phase equation. 

Assume that the initial-boundary conditions are 

(p(0,t) = (p D , E (0,1), V£G(0,T) (6.8a) 

(p(l, t) = <p Dtl (t) E (0, 1), Vt G (0, T) (6.8b) 

<p(x, 0) = <po(x), i6 0. (6.8c) 

Note here that the Dirichlet boundary condition on the right-hand side of Q is time 
dependent and the one on the left is independent of time. The problem could also 
be restated where the right-hand boundary is of Neuman type. The conditions are 



chosen to better match the experimental data that will be considered in Section 7.4 



Theorem 6.5 (Existence of Weak Solution to Diffusion Equation) Suppose that 



the following conditions hold for equations (6.6) - (6.8). 

1. Q = (0, 1) and t G (0, T) C (0, oo) 

2. ip e (0,1 - e] for all x G Q and for all t G (0, T) where < e < 1 

3. S G [e, 1 — e] and S(x) G C 1 (fi) (independent of time) where < e< 1 
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I tp D ,i e L 2 (0, T; H 1 ^)) and L°°((0, T) x Q) 



Then there exists a weak solution to (6.6) - (6.8) in the sense that tp e ipn + 
L*(0,T;HZ(n)). 



Matching equation (6.6) to the form of Alt and Luckhaus (equation (6.2a)) we 



have 



b(x,z) = (l-S(x))z, a(x,z,p)=V(x,z)p, f = 0, 



m—1. 



(6.9) 



In the conditions for Theorem 6.5 we use the parameter e to define two different sets. 
This is a small abuse of notation since ip and S need not belong to exactly the same 
set. We are simply stating that both of these functions must be bounded away from 
and 1. 



Proof: We proceed by verifying hypotheses HI - H6 of Theorem 6J2 noting the 
extension proposed in [1] to non-autonomous functions. 

HI: In 1 spatial dimension it is clear that Q is an open, bounded, and connected 
domain with Lipschitz boundary. T = {0, 1}, and H°(T) > and < T < oo. 

H2: In this case we note that b(x,z) = (1 — S(x))z. Clearly b(x, 0) = 0. Define 
$(x, z) = (1 — S(x))z 2 /2 and observe that d&/dz = (1 — S(x))z = b(x,z) and 
d 2 <S>/dz 2 = 1 - S(x) > since S(x) e (0,1). Since S G C\n) (assumption 
7^3 in the statement of the theorem) it is clear that b is measurable in the 
first component. Furthermore, b is a continuous gradient of a convex function 
in the second component. Define B(x,z) = b(x,z)z — $(x,z) + $(a;,0) = 
(l-S(x))z 2 /2. 

H3: Since a is a linear function of p it is easy to see that 



(a(x,z,p {1) )-a(x,z,p (2) )) (p (1) -p (2) ) = V{x,z) (p {1) - p {2) )' 

>C e (p^-pM) 2 



(6.10) 
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where T>(x, z) > C t for all x, z (this e-dependence reflects the choice of the 
saturation function, S(x)). Given the functional form of T> it is obvious that a 
is continuous and bounded onz6 (0, 1), p G R, and is measurable in x. Hence 
a satisfies the ellipticity condition. Simply stated, the ellipticity of the diffusion 
coefficient means that the operator in question is bounded away from zero and 
is therefore invertible. 

H4: Let z £ [e, 1 — e] and p£l. From the definition of a and /, 

\a(x,z,p)\ + |/(2)| = \a(x,z,p)\ = \V(x,z)p\ 

< Cv,e\p\ 

^(l+v^y+bi) (6.11) 

for all z, where cx> j£ is the upper bound on T>(x, z) over z. Therefore a (and /) 
satisfy the growth condition. 

H5: The left Dirichlet boundary condition is fixed in time, <p(0,t) = (fD,o- It is 
assumed that ip Dfi e L 2 (0,T; H 1 ^)) and L°°((0,T) x Q). The right Dirichlet 
boundary condition is allowed to vary in time. Assumption #4 in the statement 
of this theorem guarantees that hypothesis H5 is satisfied for this boundary 
condition. 

H6: Since b(x, z) = (1 — S(x))z it is clear that b is surjective so long as S(x) ^ 1 
and that b° = <po/(l — S(x)). That is, there exists a function y>o/(l — S(x)) 
such that b° = b(x,<p°). 

Given the final assumption in the statement of this theorem we have, in particular, 
<p D e L 1 (0,T;L oo (fi)) since L\Q) C L 2 (Vt) C L°°(fi) for sets Q of finite measure 
[38] . Therefore, from Theorem |6. 2 there exists a weak solution, cp, in the affine space 



<p D + L 2 (0, T; V) where V = {v £ H\Q) : v = on {0, 1}}. ■ 
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The uniqueness of the weak solution to (6.6) - (6.8), unfortunately, doesn't fit 



Theorem 6.3 because the diffusion operator cannot be decomposed in the manner 
required. This does not mean that the weak solution is not unique, it simply means 
that this is not the tool to prove uniqueness. This small problem is left for future 
research. 



6.2.3 Limits of the Alt and Luckhaus Theorem 

The theorem of Alt and Luckhaus does not apply to the heat transport equation 
since there are advection-type terms present in that equation that can not satisfy 



the assumed form of Theorem |6.2| The next logical direction is to see if this tool 
can be used to prove existence of the coupled saturation-humidity system at constant 
temperature. The forcing term on the right-hand side of each equation is now non- 
zero. The equations are 



as_ d_ 

dt dx 



-D{S) + Cj„)™ + C 



s) dx 



S) 



chp ^dS 
dt 



*dx 
d 



K(S)p l g 
dip 



dt dx V ' dx 






(6.12a) 
(6.12b) 



If we were to define b(z) : 



-> 



as 



b{z) 



( 



X 



i o 

-z 2 (1 - z x 



\ ( \ 



) 



vv 



one can show that there does not exist a function $ : M. 2 — > M. such that b = V$. For 
this reason we restate the equations with a consolidated form of the time derivatives 
in the second equation 



dS d 
dt dx 



-D(S) + C 



»£♦«£-*«)*)■ 


= 8?(<p,S) 


(6.13a) 


|-|:(^.S)ff) = 


= -%-(<p,S) 


(6.13b) 


U - 


= (1-S)<P 


(6.13c) 
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Solving for the relative humidity in equation (c) and substituting into equations (a) 
and (b) gives 



(i-D(S) + C> s) f x + C,l (j^) - K( S y 9 ) = er (j^, S ) 

(6.14a) 



(6.14b) 



It can be seen from this form that the coupling in the time derivatives has been 
moved to a stronger coupling with the diffusion terms. It can be shown that the as- 



sociated a(-, ■) function is not elliptic in the sense required in Theorem 6.2 Therefore 
we have determined that the Alt and Luckhaus theorems don't apply to the coupled 
system in this form. 



Equations (6.14) poses the system in a form of strong coupling known as a triangu- 
lar system. A triangular parabolic system has two equations; one parabolic equation 
with a contribution to diffusion from both dependent variables and the other with a 
contribution to diffusion from only one variable [53]. Future research into the exis- 
tence and uniqueness results will likely start here as the theory of triangular systems 
is fairly well developed and may provide a springboard to results for this problem. 

6.3 Heat Transport Equation 

In this section we consider the question of existence and uniqueness for the heat 
transport equation. This is done under the assumptions that the relative humidity 
and saturation profiles are fixed in space and time. 

If the saturation and the relative humidity are considered fixed and constant then 



the thermal transport equation (5.88c) collapses to 



dT 
P c p-fc ~ V ■ {KVT) +ph + X2 VT ■ VT = 0, (6.15) 
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where \2 is given as 

X 2 = p l c l p C l T K(S)+pW p N(S,T). 

In the absence of heat sources and if \2 is neglected we arrive at the standard heat 
equation; the existence and uniqueness results of which are well known (see any 
standard text on PDEs). It is likely that C l T ~ since, in Saito [67], the authors 
indicated that the thermal liquid flux is negligible as compared to isothermal liquid 
flux. The entropy term appearing in N, on the other hand, is likely non-negligible 
and therefore must be considered. In the case where S and tp are fixed but non- 



constant, the terms in equation (5.88c) associated with VS" and V<y? are combined 



as a source term which depends on x, t, and T. Therefore, we only need to consider 



thermal equations in the form of (6.15). If h = and Dirichlet boundary conditions 
are considered then this is the exact form of the equation considered by Rincon et 
al. in [66] : 

— -V-{a{u)Vu) + b(u)\Vu\ 2 = in Q x (0,T) (6.16a) 

u = on dQx(0,T) (6.16b) 

u(x,0)=uo(x) in Q. (6.16c) 

where we have defined u such that u ■<— T — T re f with reference temperature T re f. 
Taking h = means that we must assume that both S and <p are constant in space 
and fixed in time. This is not entirely physical, but it is a step toward a general 
existence uniqueness theory for the present equations. In this problem, a(u) is the 
diffusion coefficient, a(u) = K(u), defined either by the weighted sum of the thermal 



conductivities (equation (5.78)) or by the Johansen thermal conductivity function 
(equation (5.79)). The b function is defined as % 2 as above. 

For the Rincon existence and uniqueness theorem we consider the following hy- 
potheses: 

HI: a{u) and b{u) belong to C 1 (1R) and there are positive constants ao,ai such that 
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a < a(u) < Oi and b{u)u > 0. 



H2: There is a positive constant M > such that 



max. 



da 
du 



db 
du 



<M. 



H3: Mo € Hq(Q) fl H 2 (Q) such that ||Amo||l 2 (q) < e for some constant e > 0. 
Theorem 6.6 (Rincon et al. [66], Theorem 2.1) Under hypotheses HI - H 3 there 



exists a positive constant eo such that ifO<e<eo then the problem (6.16) admits a 
unique solution satisfying 

%. ueL 2 (0,T;H$(n)nH 2 (n)) andd t ueL 2 (0,T;H^(n)) 

u. §- V-(a(w)Vw) + 6(w)|Vw| 2 = «L 2 (Ox (0,T)) 

m. w(0) = mq 



Theorem 6.7 There exists a unique solution to equation (6.15) under the following 
conditions: 

1. u(0,t) =u(l,t) 

2. h = 

3. u(x,0) = uq G Hq(Q) n H 2 (Q) and there exists e > suc/j t/iat ||Am ||l2(q) < e 
i7ere we are using u = T — T re j for a scaled temperature (so that capital T will 



represent a finite time as in Theorem 6.6) 



Proof: We will proceed by verifying the hypotheses of Theorem 6.6 



HI: From the derivation of the heat transport equation, a(u) is a weighted sum of 
thermal conductivities from the individual phases. The particular form of the 



weighted sum comes from either equation (5.78) or (5.79), but in this scenario, 
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the saturation is presumed to be constant. Therefore, in this case a(u) is con- 
stant and is trivially C 1 (f2). The functional form of b depends on the functional 
form of the entropy and the saturation. So long as the saturation is fixed away 
from zero then b is in C^O). Furthermore, b(u) is positive so b(u)u is also 
positive for all u. 

H2: Since a is a constant, da/du = for all u. The functional form of b, on the other 
hand, is not constant so this hypothesis simply states that \2 needs to have a 
bounded first derivative. Taking the entropy term from the Darcy flux as 

" 9 = ^ ln (^) +w ' 

(from the definition of the specific heat) and defining \2 accordingly we see that 
b will have a bounded first derivative so long as u + T re f = T remains bounded 
away from 0. This is, of course, always true since T is the absolute temperature. 

H3: The third assumption of the theorem satisfies this hypothesis. 



Therefore, there exists a unique solution to (6.15) with no sources and equal Dirichlet 
boundary conditions. ■ 

6.4 Conclusion 

At this point we turn our attention to the analysis and comparisons of numerical 
solutions of the equations (both individually and coupled). The existence and unique- 
ness theory presented here is by no means complete. In particular, we are missing 
a uniqueness result for the vapor diffusion equation, the theorem used for the heat 
transport equation is very limiting with respect to boundary conditions and sources, 
and we have not mentioned results for any of the coupled systems. Many numerical 
solvers will iterate coupled systems across the equations, so an existence and unique- 
ness theory for each equation is essential to give hope that the numerical method 

converges to the solution. These results are left for future work as the ultimate crux 
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of this thesis is to justify the modeling technique against physical experimentation 
and classical models. 
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7. Numerical Analysis and Sensitivity Studies 

In this chapter we build and analyze the solution (s) to the heat and moisture 



transport model summarized in equations (5.88a) - (5.88c) with constitutive equations 



summarized in equations (5.89a) - (5.89r). To simplify matters we henceforth assume 



a 1-dimensional geometry modeling a column experiment common to soil science. 



Figure |7.1| gives a cartoon drawing of a typical column experiment with a definition 
of the geometric variable x. The grains represent a packed porous medium. Flow, 
diffusion, and heat transport are assumed to travel solely in the x direction (up or 
down) . 



9 



x 



± x = Q 



Figure 7.1: Cartoon of a 1-dimensional packed column experimental apparatus. 



In this chapter we are interested in the behavior of equations (5.88) in several 



situations related to the apparatus depicted in Figure |7.1| some physical and some 
merely hypothetical. 

1. In a drainage experiment the column is saturated with the wetting phase 
and then allowed to drain under the influence of gravity. 

Possible simplifying assumptions include: constant relative humidity and tem- 
perature. 

148 



2. In an imbibition experiment the column starts partially saturated (or dry) 
and the wetting phase is introduced either at x = 1 or x = 0. If the wetting 
phase is introduced at x = 1 then the primary force driving the liquid flow will 
be gravity, and if it is introduced at x = then the pressure head from the 
reservoir drives the flow. 

Possible simplifying assumptions include: constant relative humidity and tem- 
perature. 

3. In evaporation studies, a gradient in relative humidity is introduced between 
x = and x = 1 and relative humidity is tracked throughout the column. 
Possible simplifying assumptions include: constant temperature and/or fixed 
saturation profile. 

4. In Coupled saturation and evaporation experiments the saturation and 
relative humidity are tracked throughout the column under boundary condi- 
tions that drive both. 

Possible simplifying assumptions include: constant (or at least fixed) tempera- 
ture. 

5. In fully coupled systems we consider a heat source (typically located at x — 1) 
and boundary conditions that drive all three equations. 

In Chapter [6] we discussed the questions of existence and uniqueness of solutions 



to equations (5.88). We now turn to numerical analysis. In Sections 7.1 7.2 and 7.3 
we discuss various numerical solutions associated with the situations outlined above. 
For example, in Section 7.1[ we examine numerical solutions associated with drainage 



and imbibition experiments (types 1 and 2 above). In Section 7.4 we compare with 
a 1-dimensional column experiment outlined in Smits et al. [78J. No two- or three- 
dimensional experiments are performed in this work. 

7.1 Saturation Equation 
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In this subsection we consider the saturation equation (5.88a) with fixed and con- 



stant relative humidity and temperature and no mass transfer. That is, we consider 

These assumptions are natural in an oil-water system or simply unsaturated systems 
where the relative humidity is considered fixed experimentally. We would like to 
determine qualitative behavior of solutions to this equation under certain boundary 
conditions, experimental setups, van Genuchten parameters, and values (or functional 
forms) of r and C l s . As a first step toward this analysis let us consider dimensionless 
spatial and temporal scalings. Notice that the spatial dimension can already be 
viewed as dimensionless as seen in Figure |7.1 A characteristic time for this equation 



is t c = x c /k c = l/k c where k c = (p l gK s )//j,i is the hydraulic conductivity. Multiplying 
by t c and henceforth understanding t and x as dimensionless we get 



dS 


d / 


dt ' 


dx \ 




d 

dx 



(t„e- s 'K(S) [-p' c (S) + C< ] |?) 

(^If)-|^ s ^)- < 7 - 2 > 



In the case where r = 0, the qualitative behavior can be analyzed via the Peclet 
number; the ratio of the advective to diffusive coefficients 

?e - /J, „, - ^^ P- — -. (7.3) 



+ C l s 



-p' c (S) + C l s ( pl9i ^~ m) ) 5-( 1+1 /™) (£-l/m - 1) 

Since the diffusive coefficient depends on the dependent variable it is immediately 

clear that the Peclet number will change in time and space (in the study of linear PDEs 

the Peclet number is a fixed ratio that does not depend on the dependent variable). If 

Pe < 1 then the problem is diffusion dominated whereas if Pe > 1 then the problem is 

advection dominated. In a diffusion dominated problem we expect a smooth solution 

that spreads spatially in time, and in an advection dominated problem we expect more 

advection (transport) than smoothing. In quasilinear advection diffusion equations 
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(see a standard PDE text discussing the method of characteristics (e.g. [HH1ES]), if the 
diffusion term is not significantly weighted then the advective term may yield shock- 
type solutions. For example, if the material parameters for a particular experiment 



are located in the top right of Figure 7.2 then the diffusive term is weighted very 



small as compared to the advection and a shock is more likely to develop. That being 
said, a shock-type solution is non-physical so it is not expected in these experiments. 
This gives an indication that if a shock does occur then the parameters must be 
non-physical or the numerical method is not accurately capturing the diffusion. 

From the definition of the Peclet number it is clear that the action of C l s is to 
increase the damping of the diffusion term. Given the form of the Peclet number, 
it stands to reason that damping similar to that of C l s can be achieved by choosing 
different van Genuchten parameters. For this reason we presume for the remainder 
of this work that the effects of C l s are inseparably tied up with the effects of the 
p c — S relationship. Hence we can assume that C l s ~ 0. Recall that C l s is defined (see 



equation (5.45)) as 

^ = ^(0-^=2^ 
b dS 

and is interpreted as a wetting potential. 

With the assumption that C s ~ (or is at least inseparable experimentally from 
p' c (S)), the Peclet number becomes 

p e = ( am ) g(l+l/m)r S -l/m _ jNm _ 

\1 — rn) 

The van Genuchten parameters, a and m, are independent in this form of the Peclet 
number. Furthermore, the van Genuchten capillary pressure - saturation function is 
only one of several choices for this relationship. Other common forms are the Brooks- 
Corey and Fayer-Simmons models; each of which will have their own associated Peclet 



number. Figure |7.2| shows the nature of the Peclet number as a function of these 
parameters as well as the saturation. 
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Logio(Peclet number) for S = 0.1 



Logio(P6clct number) for S = i 




van Gcnuchtcn a 



(a) Log of Peclet number for S = 0.1 

Logi()(Peclct number) for S = 0.7 



van Genuchten a 
(b) Log of Peclet number for S = 0.4 

Logip(Peelet number) for S = 0.99 




6 8 1(1 12 11 1(1 1 

van Genuchten a 



fi 8 10 12 14 1(1 

van Genuchten a 



(c) Log of Peclet number for S = 0.7 (d) Log of Peclet number for S — 0.99 

Figure 7.2: Log of Peclet numbers for various values of saturation. The point at 
(a = 5.7, m = 0.94) indicates the values used in Smits et al. [78]. Warmer colors 
are associated with higher Peclet number and therefore associated with an advective 
solution. 



In Figure 7.2 it appears that the solutions to the saturation equation (with r = 0) 



become more diffusion dominated for smaller values of van Genuchten parameters. 
As S — > 1 the diffusion term gains more traction and hence dampens the advection. 
Of course, one cannot simply choose a set of van Genuchten parameters. Instead, the 
parameters are dictated by the material properties of the soil. In the study by Smits et 



al. [78], a = 5.7 and m ~ 0.94 (indicated by the point in Figure 7.2). In this instance, 
we expect an advection dominated solution with very little diffusive damping. This 
poses a danger numerically as it is close to the regime where shock-type solutions 
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could arise. 

The third-order term can be analyzed in a similar manner. To the author's 
knowledge there is no name for the ratio of the coefficients of the third-order term to 
the diffusive term 

„ ._ T £s re s p l gK s 



-tc P ' c (s) -p> c {s)m 

(tesk>s\ ( oim 



\ pi J \l-m 
fre S Ks\ 



o(l+l/m) / Q—l/m 



V & J 



Pe := H Pe (7.4) 



Thus the plots in Figure 7.2 are simply scaled versions of H. The question that 
remains is what effect the third-order term has on the solution. To answer this ques- 
tions we examine a few solution plots. These solutions are found using Mathematica's 
NDSolve function. This build-in command is a general differential equation solver 
handling ordinary and partial differential equations, systems of equations, vector 
equations, and stiff systems. For partial differential equations it uses a finite differ- 
ence approach to discretize the spatial variable and a version of Gear's method for 
implicit stiff time stepping following a method-of-lines approach [8 



Figure 7.3 shows a drainage experiment for various values of H = [tes^s]/ Pi- 
The initial condition is given in black. A Dirichlet boundary condition (S = So) is 
given at x = 1 and a homogeneous Neumann condition (dS/dx\ x= o = 0) is imposed 
at x = 0. Gravity points in the negative x direction, so that the liquid present in 



the column is expected to drain over time. Figure 7.3(a) shows that at earlier times 
a larger value of H gives a steeper front with plausibly physical saturation profiles. 
Non-physical, non-monotonic, results are observed for Hq = 1CT 2 as seen near x = 0.8 



in Figures 7.3(b) - 7.3(d) For values of Ho smaller than 10 2 we continue to observe 



a sharper front as compared to solutions for r = (shown in blue). 
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Drainage: a=5.7, m=0. 94118 



Drainage: a=5.7, m=0.94118 



— initial conditio 

— r=0 

— H =10- 4 

— H =10- 3 

— H„=10- 2 




(a) Saturation profiles at £ = t\ 



Drainage: a=5.7, m=0.94118 
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— r=0 

— H =10-' 

— H =10- 3 

— H„=10- 2 
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(c) Saturation profiles at £ = £3 



— initial condition 

— r=0 

— H =10- 4 
— H„=10- 3 

— H =10- 2 
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(b) Saturation profiles at t = £2 

Drainage: a=5.7, m=0.94118 



— initial condition 

— r=0 

— HfIO-" 
— H o =10- 3 

— H„=10- 2 




0.1 0.2 0.3 0.4 0.5 0.0 0.7 0.8 0.0 

X 

(d) Saturation profiles at t = £4 



Figure 7.3: Saturation profiles at various times in a drainage experiment with a 
5.7,n = 17. 



To be sure that the non-physical results observed for H = 10 -2 are not due to 

numerical noise we complete a numerical convergence test on this particular set of 

initial boundary conditions. A typical convergence test of a numerical method would 

compare against a known analytic solution, but in this case there is no known analytic 

solution. For this reason, we allow Mathematica to solve the problem using the default 

spatial and temporal tolerances and then compare solutions with fixed grids consisting 

of fewer mesh points to this solution. Mathematica's differential equation solver uses 

a finite difference approach for spatial discretization. The defaults for this scheme 

are fourth-order central differences where spatial points are on a static grid and the 
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number of grid points is chosen automatically based on the initial condition. For the 



tests shown in Figure 7.3 there were 103 grid points selected automatically. To check 
this solution, we examine the relative L 2 (0, 1) error as a function of time, 

'\S(k) — £(103) 1 1 L 2 



E (N)(t) 



\s, 



(103) || L 2 



where N is the number of spatial points. In Figure 7.4 we measure E^)(t) for N 
ranging from 20 spatial points to 100 spatial points. Notice that for any fixed small 
time the relative error decreases with increasing grid size; hence indicating numerical 
convergence at that fixed time. For dimensionless time greater than approximately 
0.05, on the other hand, the error decreases at a slower rate and there is evidence 
that the numerical method may not be converging. In all cases the relative error 



grows in time until approximately 0.25. While the bump that appears in Figure 7.3 



is certainly non-physical, Figure |7.4| seems to indicate that the numerical method is 
failing in this case and the results may not be trust-worthy for this set of parameters 
and initial boundary conditions. 
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Convergence Test for Drainage Experiment with Hq = 10 2 



10° F 




0.15 0.2 0.25 

time / t c 



Figure 7.4: Convergence test for drainage experiment depicted in Figure 7.3 N is the 



number of spatial grid points. In Figure [7T3| t 1 = 0.025£ c , t 2 = 0.050£ c , t 3 = 0.075t c 
and U = O.OlOtc 



Figure [7~5] shows an imbibition experiment for various value of Hq. As before, the 
initial condition is shown in black. For this experiment, Dirichlet boundary conditions 
are enforced at both x = and x — 1. Gravity points in the negative x direction, and 
the boundary condition at x — 1 indicates that wetting fluid is being added over time. 
For Hq = 10~ 4 and Hq = 10~ 3 we see plausibly physical results and we see sharper 
wetting fronts as in the drainage experiment. For Hq > 1CT 2 we almost immediately 
see a non-physical non-monotonicity appear at the top edge of the wetting front. 
Similar behavior was observed by Peszynska and Yi jSU] for their numerical scheme, 
and they stated 

"... we cannot speculate whether the apparent nonmonotonicity of profiles 
. . .relates to a numerical instability, or to a physical phenomenon." 



It is reasonable to ask whether this is associated with numerical noise, and Figure 7.6 



shows a convergence test similar to that shown with the drainage experiment. From 
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Figure 7.6 it appears as if the numerical method is converging under mesh refinement 
for dimensionless time approximately less than 0.1. The non-monotonicity appears 
in the region where the method should be stable so we tentatively conclude that 
this effect is not a numerical artifact. Finally, we observe that for H = 10 _1 the 
advection term has been overwhelmed by the diffusion and the third-order term and 
the numerical results are completely non-physical. 
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Figure 7.5: Saturation profiles in a imbibition experiment with a = 2.5, n — 5. 
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Convergence Test for Imbibition Experiment with Hg = 10 
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time / t c 



Figure 7.6: Convergence test for imbibition experiment depicted in Figure 7.5 In 



Figure 7.5, t x = 0.025t c , t 2 = 0.050t c , t 3 = 0.075£ c , and U = 0.010t c 



Clearly there are infinitely many choices of initial boundary conditions, and the 
results presented herein inherently depend on the conditions chosen. Similar types of 
non-physical behavior can be observed for other families of van Genuchten parameters, 
but the associated plots are excluded here for brevity. An empirical conclusion is that 
for H = (resKs)/ Hi greater than 10~ 3 possibly leads to non-physical behavior in the 
numerical solution. 

To the author's knowledge, an analysis of parameters of this type has not been 
completed in the literature. We have shown in this subsection that for reasonably 
small values of r we predict sharper fronts than with the traditional Richards' equa- 
tion. 



7.2 Vapor Diffusion Equation 

Next let us consider the vapor diffusion equation under assumptions of fixed con- 
stant temperature and a fixed saturation profile. This particular study is a bit peculiar 
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since it is unlikely that a saturation profile will remain fixed during an evaporation 
(or condensation) study. Of course, we could consider 5 = everywhere and study 
only evaporation in dry porous media, but this is also not realistic as enhancement 
models depend partly on the presence of a liquid phase. In this section we compare 



the present model proposed in Section 5.4.2 to the classical enhancement model and 
to Fickian diffusion. 

{l - S) Tt = Tx{ v ^ S) lh) (present model) (7 ' 5) 



Tt = Tx {^ s ^^ D9 Tx 



1 — S)— — — ( r]( a ){S)T{S)D g ' — ) (enhancement model) (7.6) 



(1 - S)-¥- = D 9 ^ (Fickian diffusion model). (7.7) 

Recall that rj( a )(S) is the empirical enhancement factor traditionally used, t(S) is the 
tortuosity, and D 9 is the constant Fickian vapor diffusion coefficient (see equations 



(5.71), (5.73), and obviously (7.7)) . The reader should note that we are slightly 
abusing notation given that r\ previously stood for intensive entropy and r is the 
label for the relaxation term in the saturation equation. This abuse of notation is 
contained to this section and should not cause confusion. 

Qualitatively, the shape of the diffusion curve in the x — ip plane for the present 
model is rather different than those of the enhancement and Fickian models. Fig- 



ure 



7.7 gives several snapshots of a sample diffusion experiment with enhancement 
parameter a = 25, van Genuchten parameter m = 0.9, and saturated permeability 
k,s = 1-04 x 10 _10 m 2 . Observe further that the steady state solutions are different for 
the two models. This is no surprise since the nonlinearities in the diffusion coefficient 
have different functional forms. 
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Qualitative comparison of diffusion models 



Qualitative comparison of diffusion models 
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(c) Relative humidity profiles at t = t% 
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(d) Relative humidity profiles at t = £4 



Figure 7.7: Sample diffusion experiment comparing the enhancement model to the 
present model. Here, a = 25, m = 0.9 (n = 10), and k = 10~ 10 with Dirichlet 
boundary conditions and an exponential initial profile. 



In Figure 5.5 we saw that there is potentially a marked difference between the 



diffusion coefficient in the present model and the enhancement model. Figures 7.8 



and |7.9| show a comparison of the diffusion coefficients for several values of the van 
Genuchten m parameter and two different saturated permeabilities. The functional 
dependence of the diffusion coefficient in the present model on the van Genuchten 



parameter can be readily seen between Figures 7.8(a) and 7.8(d) (similarly, 7.9(a 



and 7.9(d)), and the functional dependence on the saturated permeability can be 
seen between the two sets of figures. 
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Comparison of Diffusion Coefficients 



Comparison of Diffusion Coefficients 
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(c) Comparison for m = 0.8 (n = 5) 
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Figure 7.8: Comparison of diffusion coefficients for various van Genuchten parameters 
all taken with k s = 1.04 x 10 -10 and s = 0.334 to match the experiment in 
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Figure 7.9: Comparison of diffusion coefficients for various van Genuchten parameters 
all taken with k s = 4.0822 x 10~ n and e = 0.385 to match the experiment in [68] . 



In this thesis we propose that there is a relationship between the fitted a value 
and the material properties. 

Proposition 7.1 Given the van Genuchten m (or equivalently, n) parameter and 
the saturated permeability, k$, of the soil there is an a-priori estimate of the fitting 
parameter a. 



The immediate consequence of Propostion 7.1 is that if the fitting parameter can be 
predicted with the use of experimentation then it is, indeed, unnecessary. 
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To test Proposition 7.1 we use a simple heuristic approach to match the material 
coefficients, m and k$, to the calculated fitting parameter, a. This is done on the 
experiments by Smits et al. [75] and Sakai et al. [68]. In Smits et al., k$ = 1.04 x 10~ 10 
[m 2 ] and m = 0.944 with a statistically tuned a- value of 18.2. In Sakai et al., Ks — 
4.082 x 10~ n [m 2 ] and m ~ 0.799 with a-values of 5,8, and 15 considered. Sakai 
et al. indicated the best agreement with a = 8 while simultaneously considering a 
modified van Genuchten (Fayer-Simmons) model for the soil water retention curve. 
We do not consider the Fayer-Simmons model here, but as the Fayer-Simmons model 
is designed to give better agreement of the capillary pressure - saturation relationship 
with very low saturations we don't believe this negates our approach. 



The heuristic tests of Proposition 7.1 is as follows. The steady-state mass fluxes 



predicted by the propsed new model and the tranditional enhanced diffusion model 
are 

PsatD( m , Ks ){v,S)Vip and p sat r] {a )(S)r(S)D 9 Vp. 

Assuming that the mass fluxes are equal gives the equation 

V {m>Kg) ((p, S)Vp = 77 (a) (5)r(5)^V^. 

Making the further assumption that the gradients in relative humidity are the same 
at steady state then the diffusion coefficients must be equal. If this mass flux is taken 
at a liquid-gas interface we can assume that ip = 1. Hence, the left-hand side of this 
equation is a function of S, m, and k s while the right-hand side is a function of S 
and a 

V (m , Ks) (l,S)= V(a) (S)T(S)D9. 

At this point we could proceed by simply choosing a value for S and making com- 
parisons or we could consider the integral over all of S to remove the dependence on 
the saturation. We choose the latter as it gives a cumulative effect of the diffusion 

coefficient over the entire range of saturations. 
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Therefore, for each m and k$ and for fixed if = 1, there is a value of a such that 

(7.8) 
The left-hand side is a function of material parameters and the right-hand side is 



/ V {m>Ks) (l,S)dS= [ V(a) (S)T(S)D°dS. 
Jo Jo 



a function of a. The right-hand side of equation (7.8) integrates easily to a linear 
function of a 



V{a) (S)T(S)D°dS 



D'J f 1 

3 V 1-e 



The left-hand side of equation (7.8), on the other hand, isn't readily integrable due 



to the nonlinear nature of the van Genuchten relative permeability function. For this 



reason we seek an approximate solution to equation (7.8). 



Figure 7.10 shows the left- and right-hand sides of equation (7.8). The intersec- 



tions indicate the triple (m, Ks, a) where the equation is true, and hence indicates 
where the two models have the same cumulative diffusive effect over S G [0, 1]. For 



example, in Figure 7.10 , if m « 0.5 and k s rs 10 10 then we predict a fitting parameter 
of a « 30. 



In Figure 7.10, the blue and green curves are the right-hand sides of equation 



(7.8) for different saturated permeabilities. The blue curve is included to show the 



agreement with Smits et al. The green curve is included to show the agreement with 
Sakai et al. Observe that the experimental values are close to the values that make 



equation (7.8) true (the intersections indicated in the figure). This is to say that 



given m and Ks, equation (7.8) could have been used as an a-priori estimate of the 



value of a in these two experiments. Table 7.1 gives a more concise summary of the 



results found in Figure 7.10 
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Cumulative contribution of diffusion coefficients over S € [0, f] 
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Figure 7.10: The blue and green curves show the left-hand side of equation (7.8) for 



different saturated permeabilities, and the red lines show level curves for right-hand 
side for various values of a. The blue and green curves can be used to predict the 
value of a before experimentation. 

Table 7.1: Measured and predicted value of the fitting parameter a based on equation 



(7.8). 





a measured a predicted from (7.8) 


Smits et al. 
Sakai et al. 


18.2 w 18 
8 ^9 



Comparisons with the experiments of Smits et al. and Sakai et al. indicate that, 
while not perfect, the present model gives a diffusion equation that matches ex- 
perimental findings reasonably well without the necessity of an a-posteriori fitting 
parameter. 

7.3 Coupled Saturation and Vapor Diffusion 

In this section we couple the saturation and vapor diffusion equations under 

reasonable boundary conditions. This is done while holding the temperature fixed. 

165 



The purpose of this short study is to determine the roles of the mass transfer, ef° , the 
C dif/dx term appearing in the saturation equation, and the time rate of change of 
saturation that appears in the vapor diffusion equation. The equations are restated 
here for reference. 



(V*( S ) (-,.<< + «,|| + Ct| - A)) = J, (7.9a) 



where 

ef = M<p (p l - p sat y) ( ~ Pc + ^ S ~ pl ° - ^Tln(A^) 

This is a system of advection-diffusion-reaction equations with a pseudo-parabolic 
damping term in saturation (for r /0) and no advection term in the relative humidity 
equation. To judge the relative affect of C l compare the rate of movement of the water 
through the liquid phase with the rate of movement of water in the gas phase. As 
such, we consider the ratio of the coefficient of this term to the saturation diffusion 

Cl (Cl\ ( Clam \ 1-L1/ ,. 

~ ~r) Pe= \ * n _ S 1+1/m {S- 1/m - l) m . (7.10) 



-p'c \p l 9 \ p l g{ 1 - m ) 



In the (unlikely) case that ratio (7.10) is approximately 1 then the diffusion in relative 
humidity has equal effect as the diffusion in saturation in controlling the transient na- 
ture of the saturation. This does not fit with our physical experience so we conjecture 



that the ratio is much smaller. Figures 7.11 show time snapshots of an imbibition 



experiment with simultaneous vapor diffusion and (temperature independent) evap- 
oration. Both saturation and relative humidity are controlled with fixed Dirichlet 
boundary conditions and initial profiles consistenti with imbibition into a low satu- 
ration column. 
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Figure 7.11: Comparison of coupled saturation-diffusion models for various weights of 
Cl with parameters: k s = 1.04 x 1(T 10 , e = 0.334, H = 10" 3 , a = 4, and m = 0.667. 



From Figures 7.11, if ratio (7.10) is greater than or equal to 10 then non 



physical results are observed for this particular set of initial boundary conditions. As 
there are infmtiely many sets of initial boundary conditions we only present this one 
particular case as a proof of concept. In general we observe that this ratio must be 
kept below 10~ 3 . With a ratio this small we are simply saying that the enhancement 
in saturation seen due to increased levels of relative humidity have a small affect as 
compared to gradients in capillary pressure in the case of fixed temperature. 

7.4 Coupled Heat and Moisture Transport System 
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In this final section we examine the fully coupled system of saturation, vapor 
diffusion, and heat transfer. The culminating goal of this section is to compare 
the numerical solution of the present model with the experimental results associated 
with [78J. Dr. Smits was generous enough to share the experimental results for this 
comparison. In Section [7.4.1 the physical apparatus is discussed as well as material 



parameters and initial boundary conditions. In Section 7.4.2 the full system is solved 
numerically and compared to the experimental data. 

7.4.1 Experimental Setup, Material Parameters, and IBCs 

The experiment of interest is to track temperature, relative humidity, and satu- 
ration in a column of packed sand. Soil moisture, relative humidity, and temperature 
sensors were placed throughout a 111cm column of packed sand. A heat source was 
turned on and off above the surface of the soil (to simulate natural temperature cy- 
cles). The goal of Smits et al.was to determine whether the equilibrium assumption 
between phases was valid in porous media evaporation studies. For the our purposes 
we use this data simply as a validation of the present modeling effort. 

A schematic of the experimental apparatus used in Smits et al. [78] is shown in 



Figure 7.12 (recreated from Dr. Smits' notes). Saturation and temperature sensors 
#1 - #10, are placed every 10cm from the bottom. Saturation and temperature sensor 
#11 is 1cm under the surface of the sand. Sensor #12 is 10cm above the surface. 
Sensor #13 is on the surface ("in good contact"). Temperature sensors #14 and #15 
are placed within the insulation surrounding the apparatus (to measure the lateral 



heat loss (see the top view in Figure 7.12)). Relative humidity sensor #1 is 1cm 



under the surface, and sensor #2 is on the surface. The gray shaded area in Figure 



7.12| represents the location of the soil pack. The initial water level is the surface of 

the soil pack. The spatial variable to be used numerically is x G [0, 1] where x = 

represents the cool end of the apparatus and where x = 1 represents the surface of 

the soil 111cm above the cool end. The material properties used in this experiment 
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are shown in Table 17.21 
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Figure 7.12: Schematic of the Smits et al. experimental apparatus. Saturation and 
temperature sensors numbered 1 - 11, temperature sensors 12 - 15, and relative hu- 
midity sensors 1 and 2 [78]. The geometric x coordinate is shown on the left. (Image 
recreated with permission from 



The experiment was run for 32 days, at which point there was a power outage and 
the experiment was stopped. In the midst of the experiment there were two sensors 
that failed: saturation sensor #3 (after the 1847 tft measurement (t > 12.8days)), and 
relative humidity sensor #1 (after the 2155 ih measurement (t > 14.9days)) (see Figure 



7.13). The saturation sensors are accurate to within ±2% soil moisture content after 



soil calibration (performed by Smits et al.). The relative humidity sensor accuracy 

ranges between ±2% (for mid-range temperatures and humidities) and ±12% (for 

extreme temperatures and humidities). The temperature sensors are accurate to 

within 0.5°C for the temperature ranges of interest (www.decagon.com). 
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Table 7.2: Material parameters for experimental setup [78] . 



Parameter 


Value 


Units 


Sand Number 


30/40 


[-] 


Dry Bulk Density 


1.77 


[g cm 3 ] 


Porosity 


0.318 


["] 


Residual Water Content 


0.028 


["] 


Saturated Hydraulic Conductivity 


0.104 


[cm s x ] 


van Genuchten a 


5.7 


[m- 1 ] 


van Genuchten n (m — 1 — 1/n) 


17.8 (0.9438) 


H 



The initial and boundary conditions for the forthcoming numerical experiments 
can be taken from any point within the data set. The logical initial point for the 
numerical experiment is the beginning of the physical experiment. This particular 
point is of interest to the experimentalist as some of the interesting transient behavior 
occurs during this period. That being said, there is a significant amount of sensor 
noise in the initial phases of the experiment (see Figure 7.14(a)[ ), and if a simple proof 
of concept is all that is needed for the purposes of this work, then a later time is 
preferred so as to avoid complications related to this noise. 
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Figure 7.14: Relative humidity and temperature data showing measurement varia- 
tions in the first few days of the experiment. (Image recreated with permission from 

TO 
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Broken Sensors 
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time (days) 



Figure 7.13: Broken sensor data. Saturation sensor ^3 shown in blue and relative 
humidity sensor #1 shown in red. It is evident from these plots that these sensors 
are not working properly as they give non-physical readings. (Image recreated with 
permission from [78] ) 

The section of data where we will initially focus is between time measurements 
1800 (12.5 days) and 2150 (14.9 days). This section of data is chosen since, qual- 
itatively, it shows the least amount of sensor noise in both relative humidity and 
temperature. Saturation sensor ^3 is faulty in this time region, but the adjacent sen- 
sors indicate that there is little to no deviation from full saturation for these times. 
The relative humidity and temperature data for this time region are shown in Figure 

ma 
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Relative Humidity Data 



Temperature Data 



-C 0.5 



+3 0." 




12.6 12.8 13 13.2 13.4 13.6 13.8 14 14.2 14.4 14.6 14.8 15 



time (days) 
(a) Sensor data for relative humidity 



S-i 325 

3 



a) 

J^ 315 



-sensor #10 (11cm under surface) 
-sensor #11 (lcm under surface) 
-sensor #13 (on surface) 




12.6 12.8 13 13.2 13.4 13.6 13.8 14 14.2 14.4 14.6 14.8 15 15.2 

time (days) 
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Figure 7.15: Relative humidity and temperature data at a window beginning roughly 
12.5 days into the experiment. This window is chosen since the sensor noise is quali- 
tatively minimal in this region. (Image recreated with permission from |78j ) 



The peaks and valleys of the temperature and relative humidity data (associated 
with the on-off cycle of the heat lamp) have small variations that are likely due to 
sensor noise. To avoid modeling this noise directly we can approximate the data 
with either a simple sinusoidal function or a square wave approximation (found by 
applying the sign function to the sinusoidal approximation). The data suggests a 
square wave approximation, but the jumps in data may cause numerical difficulties 
as the derivatives at the points of discontinuity are technically delta functionals. A 



graphic of these approximations is shown in Figure 7.16 
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Relative Humidity Boundary Condition 



Temperature Boundary Condition 
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Figure 7.16: Approximations to relative humidity and temperature boundary condi- 
tions at the surface of the soil. 



Any starting point can be taken within this window of time. We choose the 
2000 th time step as the initial condition (somewhat arbitrarily) and fit functions to 
the coarse spatial data for saturation, relative humidity, and temperature. For the 
relative humidity and saturation profiles we choose hyperbolic tangent functions since 



they exhibit the primary features observed in the data (see Figures 7.17(a) and 7.17(b) 



respectively). For the temperature initial condition we choose an exponential function 



(see Figure [7.17(c) ). 



To summarize, thus far we have boundary conditions for relative humidity and 
temperature at x — 1 and we have initial conditions for all of the variables. The 
boundary conditions at x = are much simpler. For saturation and relative humidity 
we can take S(x = 0, t) = 1 and cp(x — 0, t) — 1 based on the fact that the saturation is 
fixed mechanically at 100% at the bottom end of the apparatus. For the temperature 
we can either take T(x — 0, t) — Tq or dT/dx(x = 0,t) = 0. The Dirichlet condition 
simply states that the temperature is fixed, and the Neumann condition states that 
the bottom of the apparatus is insulated so that no heat is lost. Throughout the 
course of this experiment, the thermal effects are not appreciably translated to the 
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Figure 7.17: Approximate initial conditions at the 2000 t?l data point it rj 13.9 days). 
Error bars indicate approximate sensor accuracy. 

bottom of the apparatus so either boundary condition would be sufficient. Finally, the 
only boundary condition remaining is the saturation condition at the surface of the 
soil. A simple condition is to state that the flux of liquid is zero across this boundary. 
Mathematically, this translates to the Neumann condition dS/dx(x = l,t) = 0. 

A fine point needs to be stated regarding the relative humidity equation. The sat- 
uration initial condition states that much of the experimental apparatus is completely 
saturated with liquid water (below sensor #9 approximately). The issue is that there 
is no gas phase present when S — 1. Mathematically this translates to a Stefan-type 
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problem where the lower boundary for the gas phase is actually moving spatially as 
the liquid water evaporates. For the sake of illustration let us simply assume for a 
moment that the saturation equation is a hyperbolic linear advection equation where 



the front simply advects in time. Figure 7.18 illustrates how the gas-phase domain 
might evolve in time in this simplified example. Of course, the saturation equation 
is not such a simple equation but the essence of the moving domain is the same 
regardless. 



.2 °- 6 
- 

-g 0.5 
f. 



— saturation at time 1 
— saturation at time 2 
7— saturation at time 3, 

Sil = gas phase domain at time 1 
^2 = gas phase domain at time 2 
fi 3 = gas phase domain at time 3 



0.2 0.3 0.4 0.5 0.6 0.1 

x (dimensionless length) 



n t 



Figure 7.18: Illustration of how the gas-phase domain might evolve in time. 



One way to model this Stefan problem is to assume that when S — 1 then (p — 1. 
Under this assumption we can define the gas-phase equation as a piecewise-defined 
differential equation: 



dp 



0. 



-y| ((1 - S)psat) + V ■ (p sai PVy) + ef" 



(1 - S)p. 



sat 



S = l 



S e (0,1) 



(7.11) 
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This is a somewhat artificial setup, but it allows us to assume that a relative humidity 
exists for the entire domain even when a gas phase doesn't (strictly speaking) exist. 



Simply stated, equation (7.11) indicates that if the saturation is 100% then there 
is no change in relative humidity (even though technically there is no gas phase). 
Then if we take the initial condition as ip(x, 0) = 1 when S = 1 and the left-hand 
boundary condition as ip(0, t) = 1 we (artificially) create a relationship for the relative 
humidity that holds over the entire spatial domain. There are several concerns with 
this approach, not the least of which is that the existence and uniqueness theory 
discussed previously does not cover this sort of case. Moreover, numerically requiring 
that the transition point is exactly 1 is not reasonable and some artificial cutoff, 
S* < 1, should be used to loosen this condition in numerical simulations. 

A simpler way to model the relative humidity in this experiment is to prescribe an 
initial saturation that allows for some gas phase to exist throughout the experiment 
for all time. This is achieved by setting the initial saturation less than 1. From 
a numerical standpoint this makes the equations easier to solve as the boundaries 
are all stationary in time. On the other hand, this choice of initial condition does 
not match the experimental setup and is therefore less desirable for the purposes 
of model validation. The forthcoming numerical experiments are performed using 
a combination of these two approaches; if S < S 1 * < 1, then the relative humidity 
equation is turned on and diffusion is allowed to occur. 

In the experiment, the relative humidity was only measured at the surface, at 
lcm below the surface, and in the surrounding (ambient) conditions. Hence, we can 
only truly compare with this data in regions very close to the surface. The saturation 
and temperature sensors, on the other hand, are placed coarsely throughout and the 
model can be validated over the entire spatial domain. 

7.4.2 Numerical Simulations 
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In this subsection we perform numerical simulations of the coupled system based 



on the initial and boundary conditions in Section 7.4.1 These numerical simulations 
are a first step toward validation of the newly proposed mathematical model. The 
main thrust of this work was not to create robust numerical solvers for coupled systems 
of PDEs. As such, we rely on the NDSolve package built into Mathematica as the 
primary numerical solver. The plotting and post processing are performed on a mix of 
Mathematica and MATLAB. The purpose of these numerical experiments is to provide 
a validation for the proposed equations. As such, we chose to directly model the 



Stefan problem with the relative humidity equation as shown in (7.11). 

The NDSolve package is based on a method of lines approach to numerical time 
integration with a finite difference spatial discretization. In time we use a fourth-order 
Gear's scheme and in space we use a central differencing scheme. One disadvantage 
to using this type of spatial scheme in this problem is that the saturation and heat 
equations have advective components. It is well known [54"1 ES] that upwind schemes 
are typically better at capturing the physics of advective flow problems and the central 
differencing schemes will usually introduce artificial diffusion into the solution. In 
an evaporation-type experiment such as this one, the advective flow is expected to 
be less dominant than in drainage or imbibition experiments. Hence, the artificial 
diffusion introduced with a central difference scheme is expected to have little impact 
on the solution quality. To date, the use of non-central schemes and adaptive mesh 
refinement are not supported by Mathematica's differential equation solving package. 

The parameters that can be varied in this experiment are the coefficient of the 
dynamic saturation term, r, the weight of the evaporation coefficient, M, and the 
weights of the Vy? and VT terms in the saturation equation, C and C T . There are 
no parameters that can be varied in the relative humidity equation due to the newly 
proposed model for diffusion. This is in contrast to the standard enhancement model 
where the empirical fitting parameter for diffusion is used to match the experimental 
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data. The fact that we have three different fitting parameters in the saturation 
equation simply allows us to fine tune the shape of the saturation solution beyond 



what is predicted by the traditional Richards' equation. Recall from Section 7.1 that 
larger values of the dynamic saturation term affects the sharpness of the moving front 
in the saturation equation. 



In Section 5.4.4 it was demonstrated that for certain choices of M in the evapo- 



ration rule, the present evaporation model approximated that of Bixler [19 

e^ = b(s l -e l r )R 9 "T(p sat -p 



n 9« N 



In |78j . the fitting parameter for Bixler's model was b = 2.1 x 10 5 . The corresponding 
fitted parameter in the present model is M ~ 10~ 15 where ef* is given as 

ef = M<p (ft - /?>) fa 1 - fj*>) . 

This is only an order of magnitude approximation and fine tuning can be made to 
better fit the data. 

Several experimental estimates of r were presented in Table 2 of [UJ . From this 
data, t could possibly span several orders of magnitude: 10 4 < r < 10 8 . Unfortu- 
nately, the soil types were only listed as "sand" (or "dune sand") and the relevant 
permeabilities and van Genuchten parameters were absent from this summary. These 
values at least give a ballpark estimate for experimentation with r. The coefficients of 
V<£> and VT in the saturation equation, on the other hand, are new to this study and 
appropriate values have yet to be determined. As such, we study different orders of 
magnitude for these values to estimate the effect of the terms to the overall numerical 



solution. The material parameters are all chosen to match those in Table 7.2 



The initial simulations will be run with sinusoidal boundary conditions as shown 



in Figures 7.16 This is to give a qualitative estimation of the behavior of the solu- 



tions without the trouble of the jump discontinuities associated with the square wave 

approximation or data interpolation. A smoothed square wave approximation (also 
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shown in Figures 7.16) is then used to give a closer match to the experimental data. 
The smoothing is achieved by taking piecewise defined hyperbolic tangent functions 
to approximate the steps. To measure the error between the data and the numerical 
solution we use a sum of the squares of the residual values measured at each sensor 
location: 

eu(t) = ——2 , (7.12) 

2^i££data { U data{<; , I ) ) 

where u represents any of the three dependent variables of interest (S, <p, or T) and 
the subscript indicates where the value is taken from. Stating that "£ e data" simply 
means that £ spans the sensor locations relevant for the given u (i.e. £ G {110/111, 1} 
for u = ip). Obviously e u (t) is a function of time so to get a single measure that 
describes the error we take the maximum of e u (t) over the length of an experimental 
day 

E u = max (e u (t)) . (7.13) 

A single experimental day was chosen due to numerical difficulties and due to loss of 



relative humidity sensor information. Equation (7.13) gives a single numerical value 
measuring the fit of the numerical solution to the data. In the relative humidity this 
is a very simplistic exercise as there is only 1 data point to compare against; the 
sensor located 1cm below the surface of the soil. For the saturation and temperature, 
on the other hand, this gives a better quantitative measure. 



Table 7.3 gives errors measured with equation (7.13) against the classical system 



of equations (Richards', Enhanced Diffusion, and de Vries). Table 7.4 shows the error 



as measured with equation (7.13) for various values of r, for different functional forms 



of the thermal conductivity (see Section 5.4.3), and for various values of C l and C T . 



In order to make comparisons with r, C l , and C l T we use the ratio of this coefficient as 
compared to the diffusive term in the saturation equation. This was done in Sections 
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Table 7.3: Relative errors measured using equation (7.13) for the classical mathe- 



matical model consisting of Richards' equation for saturation, the enhanced diffusion 
model for vapor diffusion, and the de Vries model for heat transport. These are com- 



pared for the two thermal conductivity functions of interest (weighted sum (5.78) and 



Cote-Konrad (5.80)). 





Relative Errors 


Conductivity 


Boundary Cond. 


Saturation 


Rel. Humidity 


Temperature 


Weighted Sum 
Cote-Konrad 


Smoothed Square 
Smoothed Square 


0.00356 
0.00357 


1.54048 
1.27818 


0.000515 
0.000631 



7.1| and |7.3[ and the ratios of interest are repeated here for clarity: 



p> c (S) 

\-Jrp 



(1^*) Pe(S) = H Pe(S) 
h|) P < S ) = R o P <S) 
(fg) P<S) = e Pe(S). 



(7.14a) 
(7.14b) 
(7.14c) 



Since each ratio is relative to the Peclet number (which is a function of S) we 



focus only on the ratios on the right-hand sides of equations (7.14). Due to the 



fact that this is a large parameter space, only some of the notable relative errors 
are presented. Mesh refinement was used in the comparisons in several instances to 
minimize numerical artifacts. Spatially, the meshes ranged between 100 and 1024 
points. Only a uniform mesh was considered. 



It is apparent in Table 7.4 that the best error approximations for saturation, rel- 
ative humidity, and temperature are found with smaller values of H (or equivalently, 
r). This observation is particular to a drainage-type experiment. If the experiment 



were an imbibition-type then it is conjectured (based on the results in Section 7.1) 



that the value of r would play a larger role. Also apparent in Table 7.4, we see that 



the values of C l and C l T play little role in the overall dynamics of the problem. 

Keep in mind that the relative humidity errors are really just the difference be- 
tween 1 single sensor and the numerical solution at that physical location. In the 

author's opinion it is unreasonable to judge the effectiveness of the numerical solu- 

180 



Table 7.4: Relative errors measured using equation (7.13) for instances within the 



parameter space consisting of the thermal conductivity function (weighted sum (5.78) 



and Cote-Konrad (5.80)), C^Cip, and r. These are taken for a (smoothed) square 
wave approximation to the boundary conditions. (The starred rows indicate failure 
of the numerical method, and the errors from the classical model are repeated for 
clarity) 



Parameters & Functions 


Relative Errors 


Conductivity 


Rq 


00 


Ho 


Saturation 


Rel. Humidity 


Temperature 


Weighted Sum 


Classical Model 


0.00356 


1.54048 


0.000515 


Weighted Sum 








IO" 25 


0.011966 


1.206005 


0.000463 








10-3.0 


0.009020 


1.201793 


0.000463 








1Q-3.5 


0.006508 


1.198274 


0.000463 








1Q-4.0 


0.004904 


1.199174 


0.000463 








1Q-4.5 


0.004076 


1.195180 


0.000463 








1Q-5.0 


0.003712 


1.196750 


0.000463 











0.003536 


1.199584 


0.000463 


Weighted Sum 


IO" 5 


IO" 5 


IO" 5 


0.003712 


1.196751 


0.000463 




1(T 4 


10" 4 


IO" 5 


0.003712 


1.196756 


0.000463 




1(T 3 


IO" 3 


io- 5 


0.003710 


1.196807 


0.000463 




io- 2 


io- 2 


io- 5 


0.003692 


1.197041 


0.000463 




io- 1 


io- 1 


io- 5 


0.003509 


1.202644 


0.000462 


-k 


1 


1 


io- 5 


* 


-k 


-k 


Cote-Konrad 


Classical Model 


0.00357 


1.27818 


0.000631 


Cote-Konrad 








1Q-2.5 


0.011964 


0.946441 


0.000516 








1Q-3.0 


0.009022 


0.951573 


0.000515 








1Q-3.5 


0.006513 


0.950024 


0.000515 








1Q-4.0 


0.004910 


0.948023 


0.000515 








1Q-4.5 


0.004084 


0.944724 


0.000516 








1Q-5.0 


0.003719 


0.946599 


0.000516 











0.003545 


0.942403 


0.000516 


Cote-Konrad 


KT 5 


IO" 5 


IO" 5 


0.003719 


0.948032 


0.000516 




10" 4 


10" 4 


IO" 5 


0.003720 


0.941801 


0.000516 




KT 3 


IO- 3 


io- 5 


0.003717 


0.945905 


0.000515 




io- 2 


io- 2 


io- 5 


0.003698 


0.947884 


0.000515 




io- 1 


io- 1 


io- 5 


0.003510 


0.966405 


0.000513 


-k 


1 


1 


io- 5 


* 


-k 


-k 



tion based solely on one point. One complication that arose within this solution is 
that the relative humidity exhibits small periods of non-physical behavior for certain 
parameters and boundary conditions. Mesh refinement removes some of this effect, 
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but even with further mesh refinement not all of the non-physical regions were re- 
moved. Possible sources of this problem are: (1) the fact that Mathematica uses cubic 
interpolation polynomials to deliver the solutions to numerical differential equations 
(cubic interpolation can overshoot sharp transitions in data), and (2) the Stefan na- 
ture of the problem causes numerical stiffness at the point of transition. Further 
studies are needed to determine the exact cause of this non-physical behavior. 



Figures 7.19 show individual time steps of several solutions for various parameters 



with a sinusoidal approximation to the relative humidity and temperature boundary 



conditions. Figures 7.20 show the same plots with smoothed square wave boundary 
conditions. The square wave boundary conditions obviously give closer approxima- 
tion to the boundary data, and at the same time the use of the square wave boundary 
conditions removes the non-physical behavior in this case. The plots associated with 
the sinusoidal approximation to the boundary conditions are presented here for com- 
parison between very smooth and slightly sharper transitions in boundary data. A 



closeup of the regions of non-physical behavior is shown in Figures 7.21(a) and 7.21(b) 
Observe in these figures that the diffusion equation solved with smaller values of H 
and larger diffusion (from the weighted sum equation) give the most plausible solu- 



tions. Figures 7.21(c) and 7.21(d)|give an indication of the difference in the relative 



humidity equations given different thermal conductivity functions. 



The comparisons of the temperature solutions are shown in Figures [7.22| and |7.23 
There is very little difference between the models for various values of r (or equiva- 
lently, Hq), so only the curves associated with Ho = 1CT 5 are shown. Observe that the 
thermal equation does a poor job capturing the extent of the diffusion near the top of 
the experimental apparatus, but it does well otherwise. Possible sources of this error 
come from: (1) the terms neglected in the simplification of the thermal model, (2) 
the initial condition, (3) the thermal conductivity functions (or parameters), and/or 
(4) the accuracy of the sensor information. 
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(d) Comparison at t = 150 x 600 sec 



Figure 7.19: Comparison of relative humidity and saturation for the fully coupled 
saturation-diffusion-temperature model as compared to data from [78]. Boundary 
conditions are taken from a sinusoidal approximation of boundary data. Thermal 



conductivities are taken as either weighted sum (5.78) or Cote-Konrad (5.80). 



Comparison of the error estimates between the two models (Table 7.3 compared 



to Table 7.4), we see that the present model gives slightly better approximations as 



measured with this metric. Table |7.5| gives the percent improvement of the present 
model over the classicl model. These values are chosen from the tables presented 
herein, and as such this is a lower bound on the percent improvement. The fact that 
there was an improvement in error is less important, in the author's opinion, than the 
fact that the models predict nearly the same error while (1) removing the necessity 
for the enhanced diffusion parameter, and (2) putting the entire system of equations 
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Figure 7.20: Comparison of relative humidity and saturation for the fully coupled 
saturation-diffusion-temperature model as compared to data from [78]. Boundary 
conditions are taken from a smoothed square wave approximation of boundary data. 



Thermal conductivities are taken as either weighted sum (5.78) or Cote-Konrad (5.80). 



on a firm thermodynamic footing. 

In this problem there are several parameters of interest, but the present study 
suggests that variations in these parameters play little role in the overall dynamics of 
the problem. This narrows us down to only 1 fitting parameter for this problem: the 
coefficient, M, in the rate of evaporation term. This was taken to best match with 
the fitted evaporation rate in [75] so it is expected that this value can be considered 
roughly constant. The classical model consisting of Richards' equation, Philip and 
de Vries diffusion equation (with enhancement fitting factors), and the de Vries heat 



184 



Saturation and Relative Humidity Comparison 
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(a) Blowup comparison at t = 100 x 600 sec. 
(sinusoidal approximated boundary) 

Saturation and Relative Humidity Comparison 



(b) Blowup comparison at t = 100 x 600 sec. 
(square wave approximated boundary) 

Saturation and Relative Humidity Comparison 




(c) Blowup comparison at t = 100 x 600 sec. 
(sinusoidal approximated boundary) 



(d) Blowup comparison at t = 100 x 600 sec. 
(square wave approximated boundary) 



Figure 7.21: Blowup comparison of relative humidity and saturation for the fully 
coupled saturation-diffusion-temperature model as compared to data from [78]. The 
inset plots give a closer look at the behavior exhibited by these particular solutions. 



transport equation contains at least two fitting parameters that are calculated using 
a least squares statistical approach. 



7.5 Conclusion 



In Sections 7.1 - 7.3, several numerical results were presented indicating the con- 
sistency of the present models with the classical mathematical models for saturation 
and relative humidity. Of particular importance is the analysis of the enhanced diffu- 



sion problem. The arguments presented in Section 7.2 indicate that modeling vapor 
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(d) Comparison at t — 150 x 600 sec 



Figure 7.22: Comparison of temperature solutions for the fully coupled saturation- 
diffusion-temperature model as compared to data from [78] . Boundary conditions are 
taken from a sinusoidal approximation of boundary data. 



diffusion in unsaturated media with the chemical potential can eliminate the neces- 
sity for a fitted enhancement factor. Also shown within this section is a sensitivity 
analysis of the r parameter for the dynamic capillary pressure term as well as the 
coefficient of Vy? that appears in the saturation equation. 



In Section 7.4 it was shown that the fully coupled system matches quantitatively 



and qualitatively to experimental data for heat and moisture transport. There are 
several problems with the matching to this experimental data. First of all, the spatial 
data is very coarse so getting an exact fit for the initial conditions is difficult. Secondly, 
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Figure 7.23: Comparison of temperature solutions for the fully coupled saturation- 
diffusion-temperature model as compared to data from [78] . Boundary conditions are 
taken from a smoothed square wave approximation of boundary data. 



the data is noisy so getting a reasonable fit for the boundary conditions (especially in 
the initial experimental times), is difficult. Lastly, the Stefan nature of this problem 
causes numerical difficulties. 

The numerical simulations presented herein indicate that the proposed models 
match both physical and experimental expectations for a heat and moisture transport 
model. The model is sensitive to the choice of thermal conductivity function and 
further investigation is needed to determine which function(s) are appropriate. The 
slight non-physical nature of the results for certain boundary conditions needs to be 
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Table 7.5: Percent improvement of the present model over the classical model usinj 



equation (7.13) as the error metric 



Boundary Cond. 


% Improvement 
Saturation Rel. Humidity Temperature 


Sinusoidal 


1.71% 22.85% 16.33% 


Sq. Wave 


1.71% 5.91% 26.78% 



investigated. Further studies (both numerical and experimental) need to be performed 
and provide a baseline for future research endeavors. The coupling of these three 
processes, all of which were derived from a thermodynamic foundation, opens to the 
door to future research endeavors on coupled processes in unsaturated soils. 



8. Conclusions and Future Work 

Throughout this thesis it has been demonstrated that HMT and the macroscale 
chemcial potential are powerful modeling tools that can be used to derive mathemat- 
ical models for rather complex phenomena in porous media. Chapter-by-chapter, the 
results are as follows. 

In Chapter [2j a short discussion of different Fickian diffusion coefficients was 
presented. While this work is not new in the sense of the creation of new theories 
or equations, it serves as an aid to understand the assumptions commonly used in 
diffusion-related research and experimentation. We showed that under certain as- 
sumptions that all of the pore-scale Fickian diffusion coefficients can be assumed 
constant. In Chapter [3j the focus turned back to porous media and the fundamental 
framework for volume averaging and Hybrid Mixture Theory were presented. This 
chapter serves as a reference for these topics and no new results were presented. 

In Chapter [4] we derived and generalized the primary constitutive equations for 
unsaturated media. In particular, we generalized the forms of Darcy's and Fick's 
laws proposed by Bennethum [12], Weinstein [81] . and others. The form of Fourier's 
law derived is similar to that of Bennethum [14], but contains terms particular to 
multiphase media. To the author's knowledge, the chemical potential form of Fourier's 
law is new to this work. The exact generalizations of Darcy's and Fick's laws presented 
here are novel to this work, but similar terms have been proposed in other works 
[12] [j~3l I8"T] . What is novel to this work is the extension of these terms to multiphase 
flow and the use of the macroscale chemical potential as a dependent variable to 
obtain additional insight. 

In Chapter |5j a coupled system of equations for heat and moisture transport 
was derived. Of particular importance are the generalizations of Richards equation 
and the Philip and de Vries vapor diffusion equation. It was demonstrated that the 
enhanced diffusion model of Philip and de Vries can be re-framed in terms of the 
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macroscale chemical potential. This re-framing removes the necessity of the enhance- 
ment factor proposed by Philip and de Vries. The coupling of the vapor diffusion 
and saturation terms is achieved through the chemical potential. The relationship 
between relative humidity and chemical potential is well known in chemistry and 
thermodynamics, but to the author's knowledge it has not been previously used in 
the porous media literature. Also in Chapter [5j a generalization of the heat transport 
equation given by Bennethum [14| to multiphase media was derived via HMT. This 
new model collapses to Bennethum's model in the case of a saturated porous medium 
and also suggests corrections to the classical de Vries model proposed in 1958. 

In Chapter [6j the questions of existence and uniqueness were studied for each 
individual equation (while holding the other dependent variables fixed). These re- 
sults were known for the saturation equation, but the results presented for the vapor 
diffusion and thermal equation are unique to this work. This is, of course, because 
these equations are novel to this work. More work needs to be done on this front. 
In particular, the uniqueness result for vapor diffusion equation is absent. Also, the 
results for the thermal equation depend on strict boundary conditions which are not 
necessarily met physically. It is emphasized that the existence and uniqueness results 
presented herein are only preliminary. 

The numerical results in Chapter [7] are performed for validation purposes. When- 
ever presenting new equations it is necessary to compare against existing models and, 



when possible, experimentally obtained data. To that end, in Sections 7.1, 7.2 and 



7.3 , numerical solutions to each equation were presented along with sensitivity analy- 



ses and comparisons to classical equations. In Section 7.4[ numerical solutions to the 



coupled system were presented and compared to experimental data. It was shown 
that under certain parameter values the newly proposed model agrees with the ex- 
perimental data. 
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There are many avenues for future research left uncompleted in this work. Re- 



laxation of the assumptions outlined at the beginning of Chapter p^ (Sections 5.2 and 



5.4) provide the initial avenues for further research. 



1. It is well known that clay soils swell when wetted. Relaxing the rigidity as- 
sumption on the solid phase would allow for this phenomenon. Mathematically, 
though, this creates further complications in the liquid- and gas-phase mass 
conservation equations since the divergence of the solid-phase velocity will no 
longer be zero due to the remaining terms in the solid-phase mass balance equa- 
tion: 

e s + e s V -V s = ^e*. 

The right-hand side of this equation may be zero in most cases (no dissolution or 
precipitation of solid particles), but a constitutive equation would be needed for 
e s or e s . This has further complications in that saturation, S = e /(l — e s ), now 
varies with solid-phase volume fraction as well as liquid phase volume fraction. 

The stress in the solid phase was discussed in Chapter [4] as a consequence of 
the entropy inequality. Upon further investigation this will give a generalization 
to the Terzhagi Stress Principle which states that the fluid phases can support 
some stress on the solid matrix. This principle will be necessary to express the 
stresses on the deforming solid. 

2. Changing the soil parameters and studying the sensitivity to empirical and 

measured relationships can be another avenue of future research. In the present 

study the soil is assumed to be isotropic. For this reason it is possible to 

take the permeability as a scalar function. This is not necessarily a reasonable 

assumption in real physical problems, and one possible avenue of research is to 

relax this assumption and take a full tensor representation of the permeability. 

Clearly a 2- or 3-dimensional simulation would have to be considered in this 
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research. 

Another adjustment to the soil properties is to the relative permeability and 
capillary pressure - saturation functions. The van Genuchten relationships were 
used for this work, but these are not the only functional forms available. The 
Fayer-Simmons model [36], for example, is an extension of the van Genuchten 
capillary pressure - saturation relationship that accounts for very small satura- 
tion content. A simple avenue for future research is to implement this model 
(and others like it) and the study to the sensitivity of the models to these small 
changes. 

3. The present model can be extended to multiple fluid phases. In this model we 
have only considered one liquid phase and one gas phase. It is reasonable to 
assume that the gas phase is a binary ideal gas, but it is not always reasonable 
to assume that the only fluid present is a pure liquid. 

One possible avenue for future research is to assume that the liquid phase has 
a dissolved contaminant. In that case it is not reasonable to assume that the 
density of the liquid is only a function of temperature. Furthermore, the dif- 
fusive term in the liquid phase equation may not be absent (depending on the 
type of contaminant). 

Another possible case is where several liquid phases are present. In this case it 
would be necessary to include a mass conservation equation (and related Darcy- 
type constitutive law) for this other fluid phase. This is a common case when 
water, oil, and gas are present within the porous matrix. This type of system 
has seen recent media attention due to the practice of fracking (driving high 
pressure water and chemicals into porous rock to release oil and natural gas). 

As mentioned previously in this chapter, an avenue for future research is to shift 

the focus from modeling to analysis. The existence and uniqueness results presented 
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in Chapter [7] are incomplete, and more research needs to be completed to give a full 
analytic description of the behavior of these equations. Of particular interest is the 
study of the fully coupled system. There are several papers discussing strongly coupled 
system of reaction diffusion equations (e.g. [31 [53]). These papers serve as a starting 
point to understanding the analysis for the coupled system. 

A third avenue for future research is to focus on the numerical method for solving 
the coupled system. The numerical solutions presented in Chapter [7] were found using 
Mathematica's NDSolve package. This is a general purpose finite difference solver 
designed to solve wide classes of ordinary and partial differential equations. Even so, 
there are problems with the methods used within that package. The foremost issue is 
the use of central differencing schemes. When solving advective (or hyperbolic-type) 
equations it is often advantageous to use upwind-type numerical schemes. This is 
not possible with the NDSolve package and causes some issues with the numerical 
solutions in advection-dominated simulations. 

Future research for the numerical method can be approached in several ways: 

1. The equations of interest in this work form a system of conservation laws, and 
as such a finite volume method is likely the best choice. Peszyhska and Yi 
[dD] derived a cell centered finite difference method and a locally conservative 
Euler-Lagrange method based on the finite difference method for the saturation 
equation with the third-order dynamical capillary pressure term. These are 
both similar to finite volume methods, and the methods derived therein can 
possibly serve as a basis for extension to the coupled system. This paper is 
of particular interest as the bulk of the numerical difficulties in the saturation 
equation arise as a result of the weight of the third-order term. 

2. There are other packages available to solve general classes of partial differential 

equations. COMSOL (Computational Multi Physics) is one such package that is 

common amongst engineering groups. This is a particularly non- mathematical 
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approach to future research, but COMSOL and other packages can be used to test 
many cases of parameters and terms suggested by the entropy inequality. 

3. The numerical simulations used to compare to the experimental data were solved 
in one spatial dimension. While this assumption is approximately valid for the 
experiment of interest, these equations need to be compared against multi- 
dimensional data. One avenue of future research is to explore the numerical 
solution to the system and compare it to a two-dimensional experiment. One 
such experiment can be found in [70] . This experiment is of interest since many 
of the parameters are the same as those used in the column experiment. 

Final Words 

The initial purpose of this work was to explore the use of the chemical potential as a 
modeling tool in porous media. It was demonstrated that the chemical potential is a 
powerful modeling tool when the underlying physical processes are diffusive in nature. 
The down sides to using the chemical potential are that it is indirectly measured and 
not widely understood. Within this work the main advantage to using the chemical 
potential was to rewrite the gas-phase diffusion equation. This allowed for the removal 
of the enhancement factor from the classical diffusion equation. 

In this work we derived three new equations that, when coupled, form a set of 
governing equations for heat and moisture transport in porous media that explains 
previously unexplained phenomena. The ideas and questions proposed within this 
chapter set up a research agenda for future years of scholarly work. 
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APPENDIX A. Microscale Nomenclature 

This appendix contains nomenclature for Part I: Pore-Scale modeling. While 
there is some overlap in notation between the two parts, this appendix allows the no- 
tation in Chapter [2] to stand alone. The equation references indicate the approximate 
first instance of the symbol (these are hyperlinked in the digital version for ease of 
use). The supporting text for these equations usually gives context and more detail. 

Superscripts, Subscripts, and Other Notations 



\Ct,. ,'t/l 



: j component of a— phase 



": a— phase 



) a,b : difference of two quantities, (■)"'' 



a: (bold symbol) vector quantity 



,: a reference quantity or a quantity evaluated at a reference state 



Latin Symbols 



c 9j : molar concentration of the j th constituent in the gas phase [mol/length 3 ] 



(2.2) 



C aj : Mass fraction of j compontnent C aj = p aj /p a [-] (2.1) 



D: diffusion coefficient [length 2 /time] (2.1) - (2.4) 



Dry-, diffusion coefficient associated with the 7-form of Fick's law [length 2 /time] 



(2.1) -(2.4) 



J 9j : flux of species j in the gas phase [mass/(length 2 -time)] (2.1) - (2.4) 
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m aj : molar mass of the j constituent in phase a [mass/mol] (2.7) 



p aj : partial pressure of species j in phase a [force/length 2 ] (2.10) 



p a : pressure of a phase [force/length ] (2.10) 



R: universal gas constant [energy/(mass-temperature)] (2.4) 



R 9j : specific gas constant for species j in the gas phase [energy/ (mass- 



temperature)] (2.3) 



T: absolute temperature (2.3) 



t: time (2.15) 



v a i\ velocity of species j in phase a [length/time] (2.1) 



v a : velocity of phase a [length/time] (2.1) 



x 9j : molar concentration of j constituent in the gas mixture [-] (2.2) 



Greek Symbols 



jj, 9v : chemical potential of water vapor in gas phase [energy /mass] (2.3) - (2.4) 



fi 9v : chemical potential of water vapor at standard temperature and pressure 



[energy /mass] (2.7) 



p° J : mass density of species j in phase a [mass/length 3 ] (2.1) 



p a : mass density of phase a [mass/length 3 ] ( |2.1 ) 



p S at'- mass density of water vapor under saturated conditions [mass/length 3 



(paragraph before (2.9)) 
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APPENDIX B. Macroscale Appendix 
B.l Nomenclature 

This appendix contains nomenclature for Part II: Macroscale modeling. While 
there is some overlap in notation between the two parts, this appendix allows Part II 
to stand alone. This appendix also clarifies any notational discrepancies between the 
pore-scale and macroscale models. The equation references indicate the approximate 
first instance of the symbol (these are hyperlinked in the digital version for ease of 
use). The supporting text for these equations usually gives context and more detail. 

Superscripts, Subscripts, and Other Notations 

(■) aj : j th component of a— phase on macroscale 
(■) a : a— phase on macroscale 

(•): denotes exchange from other interface, phase, or component 
(■) a ' b : difference of two quantities, (■) a ' b = (-) a — (-) b 
(•)■?: pore scale property of component j 
a: (bold symbol) vector quantity 
A: second order tensor (matrix) 
(■)o or (■)*: reference state 
Latin Symbols 



a: fitting parameter for enhanced diffusion model [-] (5.71) 



b a \b a : External entropy source [energy/ (mass-time-temperature)] (3.37) 



197 



C a] : Mass fraction of j compontnent [-] (3.23) 



C u : Coefficients of V« terms in generalized Darcy's law (u = S, T, ip) (5.44) 



C_ s : Right Cauchy-Green tensor of the solid phase = (F S ) T ■ F_ s [-] (4.4) 



C : Modified right Cauchy-Green tensor = (F) T ■ F [-] (4.4) 



D a : diffusivity tensor [length 2 /time] (4.81) 



T>: generalized diffusivity function [length /time] (5.69a) 



d a : Rate of deformation tensor = (Vv a ) sym [1/time] (3.40) 



e aj } e a : energy density [energy/mass] (3.34) 



eJ : rate of mass transfer from phase j3 to component j in phase a per unit 



mass density [1/time] (3.21) 



e%: rate of mass transfer from phase /3 to phase a per unit mass density 



[1/time] (3.24) 



F_ s : Deformation gradient of the solid phase [-] (4.5) 



F_ : Modified deformation gradient of the solid phase [-] (4.5) 



g, g a \g a : gravity [length/time 2 ] (3.14) 



h a \h a \ external supply of energy [energy / (mass-time)] (3.34) 



i J : rate of momentum gain due to interaction with other species within the 



same phase per unit mass density [force/mass] (3.28) 



J s : Jacobian of the solid phase [-] (4.3) 



K a : hydraulic conductivity tensor for a phase [length/time] (4.36) 
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K : thermal conductivity tensor [energy/(mass-time-temperature)] (4.39) 



m: van Genuchten parameter (m = 1 — 1/n) [-] (5.51a) 



n: van Genuchten pore-size distribution parameter (n = 1/(1 — m)) [-] (5.51a) 



p a : classical pressure in the a phase [force /length ] (4.29) 



p a (fi) : cross coupling classical pressure [force/length 2 ] (4.41) 



Pc — P 9 — p ■ capillary pressure [force /length ] (5.53) 



p a : thermodynamic pressure [force /length ] (4.50) 



p a (P); cross coupling thermodynamic pressure [force/length 2 ] (4.47) 



q a : heat flux for a phase [energy/(length 2 -time)] (3.35) 



q: total heat flux [energy/ (length 2 -time)] (|4.87f) 



q a : Darcy flux for a phase [length/time] (4.67) 



Q aj : rate of energy gain due to interaction with other species within the 
same phase per unit mass density not due to mass or momentum transfer 



[energy/ (mass-time)] (3.34) 



Qa'- energy transfer from phase a to constituent j in phase (3 per unit mass 



density not due to mass or momentum transfer [energy/ (mass-time)] (3.34) 



r: microscale spatial variable [length] (|3.2) 



f Qj : rate of mass gain due to interaction with other species within the same 



phase per unit mass density [1/time] (3.21) 



R: gas constant [energy/ (mol-temperature)] (5.48) 
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R 9j : specific gas constant for j constituent in gas phase [energy/ (mass 



temperature)] (4.82) 



R a : resistivity tensor [time/length] ( |4.36 ) 



S = e /(e): liquid saturation [-] (5.2) 



t: time 



T: absolute temperature (3.39) 



t aj : partial stress tensor for the j constituent in the a phase [force/length 



(3.28) 



t a : total stress tensor for a phase [force/length 2 ] (3.29) 



To : rate of momentum transfer through mechanical interactions from phase 



(3 to the j th constituent of phase a [force/length 3 ] (3.28) 



T B : rate of momentum transfer through mechanical interactions from phase 



(3 to phase a [force/length 3 ] (3.29) 



v a : velocity of the a phase relative to a fixed coordinate system [length/time] 



(3.24) 



v aj,a _ v a.j _ ^a. diffusive velocity [length/time] (3.40) 



v a,s _ v a _ v s. velocity relative to solid phase velocity [length/time] (3.40) 



w a/ 3 . : velocity of constituent j at interface between phases a and (3 [length/time] 



(3.10) 



Greek Symbols 



a: phase 
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a: van Genuchten parameter [-] (5.61) 



j3: phase 



S(t): Dirac delta function [-] 



5A a a\ Portion of the a/3— interface in REV (3.10) 



e a : volumetric content of a phase per volume of REV [-] (3.15) 



s = e + e 9 : porosity [-] (4.1) 



7] a : specific entropy of a phase [energy / (mass - temperature)] (3.38) 



f) aj : entropy gain due to interaction with other species within the same phase 



per unit mass density [energy/(mass-time-temperature)] (3.38) 



77: enhancement factor for diffusion [-] (5.70) 



T a : macroscale Gibbs potential [energy /mass] (4.61) 



7°: indicator function which is 1 if in phase a and zero otherwise (3.10) 



k: permeability [length ] (4.68) 



A a : rate of entropy production for a phase [entropy / time] (3.38) 



X aj : Lagrange multiplier for the continuity equation of phase a (4.12) 



A Qjv : Lagrange multiplier for the N term dependence relation of the compo- 



nents in phase a (4.12) 



A = p 9v /p 9 : ratio of partial pressure to bulk pressure in gas phase [-] (5.43) 



jx a : dynamic viscosity for phase a [force-time] ( |4.68 ) 



/i Qj : macroscale chemical potential of j th species in a phase [energy /mass] 



d456| 
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u a : kinematic viscosity for phase a (4.69) 



$ a : Entropy transfer through mechanical interactions from phase /3 to phase 



a per unit mass [energy/ (mass-time-temperature)] (3.38) 



7r Q : wetting potential of a phase [force/length 2 ] (4.51) 



n a w : cross coupling wetting potential [force / length 2 ] (4.48) 



ip: relative humidity [-] (5.40) 



ip a : specific Helmholtz potential of the a phase [energy /mass] (3.39) 



p a : mass density of a phase (mass a per volume of a) [mass / length 3 ] (3.22) 



p aj : mass density of j constituent in a phase (mass «j per volume a) [mass 



/length 3 ] (3.22) 



p sat : saturated vapor density [force / length ] (5.41) 



r: tortuosity [-] (5.70) 



r: scaling coefficient for dynamic saturation term [-] (4.53) 



B.2 Upscaled Definitions 

Definitions of bulk phase, species, and averaged variables resulting from upscal- 
ing. Recall that an overbar indicates a mass averaged quantity and angular brackets 
indicate a volume averaged quantity. 

1 f , . . . ^a 1 



((■)) C 



\SVr. 



{■)j a dv and 



a\ JSV 



{p) a \SV a \Jsv c 



■)ladv 



b aj = W 



N 



b a = J2c a >b c 



(B.l) 

(B.2) 
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C a i 



e a > 



ra 1 : ra 1 r. ■ n- 

e 3 + - V 3 . V 3 V a J . V a J 

2 2 



c /3 

9 a 

h" 



N 
3=1 

e a 
\6V*\ 

N 

E# 

i=i 

— rQ 

AT 

3=1 

—a 



a o I e a i _|_ _ v a j. Q . w "j." 



/ p J («? Q/3j - v 3 ) ■ n a da 

JA al3 



g 



W + gi ■ v 3 - g a] ■ v a > 



N 



3=1 



e"p" 3 [i + r 3 v 3 -v^r 3 



(q J ) a + (t j ■ v 3 ) a - fi ■ v aj + p a *v aj (e aj + -v a > ■ v aj j 



//'■?>' | e 3 + -v 3 ■ v 3 ) . 



AT 

E 

3=1 



q <*j + t a > ■ V a ^ a - p a J ( e a > + -v aj ' a ■ v aj ' a ) v a ^ n 



Qa 3 = £ a p aj 



Qj _|_ f . v j _ U + fj v j _ v <*i?j . v 



a / — ret 



Q\ 



+r 3 ( e 3 + - V j . V 3 \ - r i e-? + -iP • v 3 



1 , 



|<*K| 



a/3 



qr' + t J • v 3 + p J I e J H — v 7 • t;- 7 J (ty a ^ — v 7 ') 



(B.3) 
(B.4) 

(B.5) 

(B.6) 

(B.7) 

(B.8) 
(B.9) 

(B.10) 

(b.ii; 

(B.12) 



(B.13) 
(B.14) 



(B.15) 



n a da. 
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J A afi J 



(B.16) 
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3=1 
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N 
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N 

3=1 
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N 
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(B.17) 

(B.18) 
(B.19) 

^B.20) 



(B.21) 

^B.22) 

^B.23) 

^B.24) 
^B.25) 

^B.26) 

^B.27) 

^B.28) 

^B.29) 

^B.30) 
(B.31) 

^B.32) 
^B.33) 



B.3 Identities Needed to Obtain Inquality 3.40 
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(B.34) 
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(B.36) 

(B.37) 

(B.38) 



(B.39) 
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APPENDIX C. Exploitation of the Entropy Inequality - An Abstract 
Perspective 

This short appendix is meant to give a brief and abstract description of how the 
exploitation of the entropy inquality works. By "abstract" we mean that we will not 
assign any physical meaning to the variables. Instead we will simply state how the 
variables relate to each other and how they relate to the full set of chosen independent 
variables. The secondary purpose of this appendix is to make clear a few assumptions 
related to constitutive equations that are necessary in order for the exploitation of 
the entropy inequality to be successfull. We conclude with an inequality that dictates 
how the linearization of the constitutive relations must behave in order not to violate 
the second law of thermodynamics. This is similar to the 1968 Nobel Prize winning 
analysis by Onsager, who showed the reciprocal relations that must hold at equilibrium 
for irreversible processes. 

Let S be the set of all independent variables for the Helmholtz potential. This 
set defines the physical system of interest. Define the following sets: 

• i x j} := the set of all variables that are neither constitutive nor independent. 
Examples typically include T, VT, Vp^, 

• {Vk} '■= the set of all constitutive variables which are zero at equilibrium. Ex- 
amples typically include e a , e a J , and r a K 

• {Vk,} '■= the set of all constitutive variables which are not zero at equilibrium. 

(Xj 





Examples typically include T J , and i 3 . 



{z{\ := the set of all variables that are zero at equilibrium. Examples typically 
include VT, v a ' s , v a ^ a , and d°. 

206 



Since each z\ is an independent variable it is clear that {z{\ C S. Furthermore we 
observe that {xj} fl S = 0. The constitutive variables, on the other hand, are known 
to be functions of variables in S and as such the statement that ({yk} U {y K }) fl S = 
is easily misinterpreted. It is a true statement that ({yk} U {y K }) fl 5 = 0, and it is 
correct not to choose constitutive variables as independent variables. The confusion 
is in the fact that y k = Vk{S) for all k (and for all k). 

The rate of entropy generation, A, can be written as a linear combination of the 
variables from {xj}, {yk}, and {zi}, where the coefficients are functions of variables 
from S. That is, 

o < J2 K ° = A = E ( x 3 x i) + E (y^ + E (^ z o ^ ° ( c - x ) 

a j A; I 

where 

X d = Xj(S,Y K (S)), Y k = Y k (S,Y K (S)), and Z, = ^(^f^)). (C.2) 

This is not the only way to algebraically rearrange A, but this is what is commonly 
done during the exploitation process. 



We now use inequality (C.l) to derive equations that hold for all time, at equi- 



librium, and near equilibrium. 

C.l Results that Hold For All Time 

We have no control over the variables Xj since they are neither constitutive nor 
independent. This means that they could be positive or negative, large or small. In 
order for the second law of thermodynamics to hold for all time, the coefficients Xj 
must therefore be zero for all time. This implies that 

k I 

C.2 Equilibrium Results 
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The definition of equilibrium is the state at which all of the variables {y k } are 
zero. This definition is based on physical intuition and will vary depending on the 
system of interest. From thermodynamics, the rate of entropy generation must be 
minimized at equilibrium. This implies that the gradient of the entropy generation 
function must be the zero vector (as understood with S as the independent variables 
for the gradient). 

<9A 







dsi 



for all i. 



(C.4) 



ry 



Taking this partial derivative of the right-hand side of (|C3|) we see that 
dk 







dsi 



eq 



?(»£) + ?(* 



dyt 
dsi 



E 



Zl 



dsi 



E h 



dzi 
dsi 



eq 

(C.5) 



Since zi\ eq = = yk\eq for all I, k and ^ = 5u (since zi G S VZ) we get 







dA 
ds. 



"i 



Eh 



dy k 
' dsi 



ZiS. 



iOil 



(C.6) 



<<i 



At this point we make an assumption that greatly affects the constitutive vari- 
ables. 



Assumption: At equilibrium we must have 



dyk 
dzi 



for all L k 



(C.7) 



eq 



Under this assumption it is clear that Z\ = at equilibrium for all I. Notice that 
this says nothing about when Sj ^ {z{\. From this argument, each equation Z\ = 
gives a constraint on some of the variables in S. The assumption made can be viewed 
as a further restriction on the constitutive variables, but it is not clear whether this 
assumption is physical. 



C.3 Near Equilibrium Results 
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For the near equilibrium results we consider two types of variables: variables that 
are zero at equilibrium and constitutive variables. In each case we linearize about 
the equilibrium state. A typical linearization result for variables which are zero at 
equilibrium is 

dZ t 



\n.eq 



\Zi)\ eq -\- y j — z m 

■ ^- s — OZ m I eq 

=0 



+ £ 



dy„ 



Vn + 



"i 



s^dzn ^ dz l I 

OZ m I eq Oy n I eq 



(C. 



The value of [Z\)\ eq is zero by the above arguments, and the partial derivatives are 
now functions of all of the other variables which are not zero at equilibrium: 



D ln :-- 



dZ t 



oz m 
dZ t 



dy v 



eq dz m 

eq dy n 



(6,6,- •• 



(■(/ 



(6,6, ••■ 



(C.9) 
(CIO) 



'■■<! 



for £„ G S \ {zi}. Written more simply 



Zi\ 



n.cq 



^ImZm ~l~ i-J\nyn 



(c.11; 



where the summations are implicit over repeated indices. 

For the constitutive variables we do a similar linearization, but note that the 
equilibrium state is not necessarily zero. Therefore, 



dY k 



Y k \n.eq- (Y K )\ eq + 2_^—— y p + >J 
„ Oy P eq I —f 



dZn 



Zn + 



<<1 



Making similar definitions as before, 



Ek P ■ 



1 kq 



dY k 



dy P 
dY k 



dz n 



<<i 



"i 



the linearization result for the constitutive equations is 



(C.12) 



(C.13) 
(C.14) 



Y k 



n.eq 



1 k\eq ~r -C-'fcpi/p ~T~ "kqZq 



(C.15) 
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where the summations are implicit over repeated indices. 

The trouble here is that we must have some information about the equilibrium 
state of the constitutive variable. The presumption that this is zero may be non- 
physical. An example of this is the capillary pressure relationship derived in multi- 
phase unsaturated media. 



Pc = (Pc)\e q + re 



(C.16) 



where p c is the capillary pressure, r = dp c /di , and the equilibrium capillary pressure 
is given as a function of saturation p c \ eq = p c {S) via the van Genuchten approximation. 
This equilibrium constitutive equation is known not to be zero. In other systems 
the issue may be more subtle, but in any case one needs to have some information 
(whether from experiments or from other theory) to define the equilibrium state of 
the constitutive variable. 



C.4 Linearization and Entropy 



Consider again equation (C.3), but now substitute the linearized results into Y k 
and Zi 

< A = (y k Y k ) + { Zl Z x ) 

= Vk (Y k \ eq + E kp y p + F kq z q ) + zi (Q m z m + D hn y n ) 

= VkYk\e.q + VkEkpVp + VkFkqZq + Z/C/ m Z m + Z\Di n y n (C.17) 

(summations are again taken over repeated indices). Recognizing the quadratic terms 
as matrix products and rewriting in block matrix form gives 

\ 



<ykY k \eq + 



'■y<BF\ 



w 



\m 



ykY k \ eq + C T AC 



(C.18) 



Notice that in the absence of constitutive variables (C.18) simplifies to 



< z 1 Cz. 



(C.19) 
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Simply stated this means that G must be positive semidefinite in order for the second 
law of thermodynamics to hold. This is Onsager's Nobel Prize winning result. Recall 
that Zi\ neq = C\ m z m . If we take, for example, z m = VT and observe that Zi\ neq 
is minus the heat flux near equilibrium (the coefficient in the entropy inequality 
associated with VT is minus the heat flux), then the / — m entry in G is the heat 
flux tensor. Onsager's result dictates the positivity of the heat flux tensor and the 
linearized result give Fourier's law: —q = .K~VT. In other words, there is no accident 
that many physical "laws" take the same form as Fourier's law; they are a result of 
the non-negativity of G and the entropy inequality. This is a simple example, but 
it should help to elucidate the problem that arises when constitutive equations are 
introduced. 



Returning to equation (C.18) we see that it is not immediately obvious that A. 
needs to be positive semidefinite. In fact, the only way that we can guarantee that 
A. has this property is if ykYk\ eq < 0. That is, there must be a physical restriction on 
VkYk\eq that, when violated, one perceives nonphysical results. 

To give a physical example of this we return to the capillary pressure example. 
In this case we have y^ = i and Yk\ eq = p c {S) where we are ignoring all other 
constitutive equations (or we are taking Y r \ eq = for all r ^ k). The restriction 
derived herein states that p c (S)e l < for all time. The time derivative can clearly 
take either sign, but what this seems to be indicating is that in drainage (when e l < 0) 
the equilibrium capillary pressure must be positive, and in imbibition (when e l > 0) 
the equilibrium capillary pressure must be negative. This is a bit contradictory since 
"drainage" and "imbibition" are non-equilibrium phenomena, and as such it is not 
possible to measure the equilibrium capillary pressure at these states. This leaves us 
with a conundrum: Is there a fundamental misinterpretation of the capillary pressure 
in this example, or is there is a deep-seated flaw in the exploitation of the entropy 
inequality. 
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APPENDIX D. Summary of Entropy Inequality Results 

The following is a concise collection of the results derived from the entropy in- 
equality. This appendix is to be used for reference when building the macroscale 
models. 



D.l Results that Hold For All Time 



Helmholtz potential and entropy are conjugate variables (equation (4.14)) 

dip c 



dT 



-v 



(d.i; 



Lagrange multiplier for fluid phase (equation (4.15)) 



A* = j: 






(D.2) 



Lagrange multiplier for the dependence of the diffusive velocities on the N 



tii 



species (equation (4.16)) 



1 N 



J=l 



Solid phase pressure (equation (4.19) 



P° 



^(r) = -^£Uv^ 






dJ s 



(D.3) 



(D.4) 



Solid phase stress (equation (4.23)) 



t s = -p s I + t s + -t[ + -t{ 



:, £ s= h £ s= h 
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(D.6a) 
212 



~- h \ H = dc 
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(D.6c) 



D.2 Equilibrium Results 



Fluid pressures (equations (4.29), (4.41), (4.42), and (4.47) - (4.52)) 
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Momemtum transfer between phases (equation (4.31)) 
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dJ s dC v= ' 



(D.16) 



Momentum transfer between constiuents (equation (4.34)) 



Y^T a / + T 3 = -V (^f) + A Q ^ V (e a p a i) + V a V (5 Q p Qj ) . (D.17) 



Partial heat flux (equation (4.35)) 



J2e a q a = 



(D.18) 



Chemical potential definition (equations (4.56) and (4.57)) 
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(D.20) 



Mass transfer (equation (4.60)) 
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D.3 Near Equilibrium Results 



Momentum transfer between phases (equation (4.36)) 



E* 

va^/3 / near W/3 



W/* / eq 



(D.22) 



Momentum transfer between constituents (equation (4.38)) 



E* 



a j , ^j 



felt + H 

We* / 



- e a p a iR a i ■ v a i' a . (D.23) 
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Partial heat flux (equation (4.39)) 



E £< V | -K-VT 



(D.24) 
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D.4 Constitutive Equations 

• Darcy's law 



Pressure formulation (equation (4.66)) 



e^g ■ (VV' S ) 



N 



+ £ 



^'^-V^)v/. 






eV t^ VJS - eV lf : v s S) + v • {i ■ £) ■ (D - 25) 



Chemical potential formulation (equation (4.76)) 



H-(eV s ) 



N 



J2 (/' V/0 - pVVT + A + lv • (V : ^] (D.26) 



i=i 



Fick's law (equation ( |4.80| )) 



e a p a >R a i ■ v a >' a = -p^Vfi * + p a ig. 



(D.27) 



Total heat flux (equation (4.92)) 



a=l,g 



N 






(D.28) 
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APPENDIX E. Dimensional Quantities 

This appendix contains typical values for the quantities found in the macroscale 
heat and moisture transport model. Note that since many of the tables are large so 
some are turned sideways and some are bumped to different pages by default. 

Table E.l: Dimensional quantities 



Symbol 


Quantity 


Dimensions 


reference value 








liquid (water) 


gas (air) 


e a 


volume fraction 


— 


— 


— 


P a 


density 


ML~ 3 


1000^ 

mr 


1 k 9 


p 9v 


vapor density 


ML' 3 




2 x 1(T 2 ^ 

mr 


p 9a 


dry air density 


ML' 3 


— 


1.2-*% 

mr 


T 


temperature 


K 


298.15if 


298.15if 


B?° 


gas constant (air) 


L 2 r 2 K- 1 


— 


286-9^ 


R 9v 


gas constant (vapor) 


LH- 2 K- 1 


— 


461-5^ 


D a 


diffusion coefficient 


LH- 1 


2 x 1(T 9 ^ 


2.5 x 1(T 5 ^ 


fl 9 " 


chem. potential (vapor) 


LH- 2 




-1.27 x 10 7 f 

kg 


H 9a 


chem. potential (air) 


L 2 t- 2 


— 


1.47 x 10 7 f 

kg 


9 


gravity 


Lt- 2 




9.81*1 


K 


permeability 


L 2 


(see Tab. 


E.2 


) 


(see Tab. 


E.2) 


Ha 


dynamic viscosity 


ML- l r l 


lO-^Pa ■ s 


lO~ b Pa ■ s 


rf 


specific entropy 


L 2 r 2 K- 1 


3886 kg J K 


6519 . J „ 

kg-K 


M 


evaporation coefficient 


tL~ 2 
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